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Zusammenfassung 



Diese Doktorarbeit behandelt die Modellierung und Simulation von Heteroepitaxie- 
Wachstum. Dabei wird insbesondere der Gitterunterschied der am Wachstumsprozess 
beteiligten Materialien miteinbezogen. 

Im Einleitungskapitel wird ein Uberblick iiber die wichtigsten Oberflachenprozesse 
und die wesentlichen Mechanismen des Verspannungsabbaus beim heteroepitaktischen 
Wachstum gegeben. Es folgt eine Zusammenstellung gangiger Methoden der Modellierung 
und Simulation von Heteroepitaxie wie sie z.B. in der Molekularstrahlepitaxie realisiert 
ist. 

In Kapitel 2 wird die so genannte Molekular-Statik Methode zur Berechnung von 
Diffusionsbarrieren vorgestellt. Unter der Annahme, dass die Teilchen des Kristalls iiber 
Paarpotentiale miteinander wechselwirken berechnen wir Hupf- und Austauschbarrieren 
fur Diffusionsereignisse am Rand von Adsorbatinseln. Hierbei interessiert vor allem 
die Abhangigkeit der Barrieren von der Inselgrofie, vom Gitterunterschied und dem 
verwendeten Potential. Es zeigt sich, dass die Barriere fur Austauschdiffusion besonders 
stark von diesen Faktoren beeinfTusst wird. Im Bereich groBer Gitterunterschiede zum 
Beispiel nimmt diese Barriere deutlich ab. Damit wird hier die Austauschdiffusion zum 
dominanten Mechanismus der Interlagendiffusion. 

Im folgenden Kapitel 3 wird ein von uns weiterentwickelter Algorithmus zur Simulation 
von Heteroepitaxie- Wachstum vorgestellt. Um die Verspannungen und ihren Abbau 
realistisch modellieren zu konnen, bedarf die Beschreibung von Heteroepitaxie der 
Moglichkeit kontinuierlicher Abstande zwischen den Teilchen des Kristalls. Der Leit- 
gedanke dieser gitterfreien Methode ist es daher, die Diffusionsbarrieren fiir jedes Teilchen 
auf der Kristalloberflache unter Verwendung von Paarpotentialen als Funktion der 
Teilchenabstande zu berechnen. Die so gewonnenen Aktivierungsenergien werden dann 
in einer Kinetischen Monte Carlo (KMC) Simulation zur Modellierung der Wachstumsdynamik 
verwendet. Als entscheidender Vorteil von KMC gegeniiber anderen Methoden (wie z.B. 
der Molekular-Dynamik) konnen so die relevanten Zeitskalen des Kristallwachstums 
abgedeckt werden. Wir beschreiben in diesem Kapitel ausfuhrlich die Berechnung der 
Diffusionsbarrieren, die Relaxation des Kristalls sowie Strategien fiir eine effiziente 
Umsetzung des Algorithmus. 

Die folgenden Kapitel beschaftigen sich mit der Analyse von drei wichtigen Mechanismen 
des Verspannungsabbaus unter Verwendung der gitterfreien KMC Methode. In Kapitel 
4 untersuchen wir den Verspannungsabbau aufgrund der Bildung von Versetzungen im 
Adsorbatfilm. Im ersten Teil des Kapitels werden Mechanismen der Versetzungsbildung 
ausfuhrlich anhand von Simulationsergebnissen diskutiert. Weiterhin wird die Abhangigkeit 
der so genannten kritischen Filmdicke fiir das Auftreten von Versetzungen in Abhangigkeit 
vom Gitterunterschied untersucht. Wir finden dabei in Ubereinstimmung mit zahlreichen 
experimentellen Ergebnissen, dass sich die kritische Filmdicke mit Hilfe eines Potenzgesetzes 
als Funktion des Gitterunterschieds beschreiben lasst. Der zweite Teil des Kapitels 
behandelt das pseudomorphe Wachstum mit anschliefiender gradueller Relaxation des 
Adsorbatfilms fiir den Bereich relativ kleiner Gitterunterschiede. Diese Untersuchung ist 



durch neuartige in-situ Messungen der vertikalen Gitterkonstante am System ZnSe/GaAs 
motiviert. Es zeigt sich eine sehr gute qualitative Ubereinstimmung der Simulationsergebnisse 
mit dem Experiment. AuBerdem kann gezeigt werden, dass der Bereich pseudomorphen 
Wachstums naherungsweise nach einem Potenzgesetz mit dem Gitterunterschied skaliert. 

Im folgenden Kapitel 5 gehen wir auf die Entstehung selbstbildenderliaseln (Stranski- 
Krastanov Wachstum) als einen weiteren moglichen Relaxationsmechanismus der Heteroepitaxie 
ein. Wir fuhren hier die Bildung einer benetzenden Adsorbatschicht auf verlangsamte 
Adsorbat-Diffusion auf dem Substrat zuriick. Die anschliefiende Inselbildung finden 
wir durch zwei kinetische Faktoren begiinstigt: eine verlangsamte Diffusion auf der 
Inseloberflache und eine gerichtete Diffusion zur Inselmitte hin. Beide Phanomene haben 
ihren Ursprung in der teilweisen Relaxation von Adsorbatmaterial in den Inseloberflachen. 
Weiter bestimmen wir die Abhangigkeit der InselgroBe und Inseldichte vom Gitterunterschied, 
der Temperatur und dem Teilchenfluss. Dabei werden sehr gute qualitative Ubereinstimmungen 
mit MOVPE (metal-organic vapor phase epitaxy) Experimenten erzielt. 

1m abschlieBenden Kapitel 6 untersuchen wir anhand eines Dreikomponentensystems 
die Bildung von Oberflachen-Legzerungen als moglichen Relaxationsmechanismus. Anhand 
von Gleichgewichtssimulationen gelingt es nachzuweisen, dass die Konkurrenz von Teilchenbindungen 
und Gitterunterschied zu einem regelmaBigen Streifenmuster der beteiligten Materialien 
fiihrt. Wir untersuchen anschlieBend inwieweit sich diese Musterbildung in ein kinetisches 
Modell iibertragen lasst. Abhangig von Temperatur, Gitterunterschied und Potentialtyp 
finden wir sowohl die Streifenbildung als auch die experimentell berichtete Verastelung 
der Oberlachenstrukturen. Der abschlieBende Vergleich mit einem Gittergasmodell zeigt, 
dass Musterbildung zwar allein aufgrund kinetischer Ursachen moglich ist, der Gitterunterschied 
zwischen den beteiligten Materialien aber die Verastelung und Stabilisierung der Strukturen 
bewirkt. 

SchlieBlich werden im Anhang der Arbeit Ursprung und Eigenschaften der verwendeten 
Potentiale fur die Teilchenwechselwirkungen besprochen. 
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Abstract 



In this PhD thesis,we develop models of heteroepitaxial growth, where the lattice misfit 
of the involved particle species is of special interest. In the introductory chapter 1 we 
introduce important physical processes which occur on the crystal surface. We give an 
overview on relevant strain relaxation mechanisms and discuss different methods for 
the simulation of heteroepitaxial growth. 

In chapter 2, we introduce the so-called Molecular Static method for the calculation 
of diffusion barriers. Provided that the particles of the crystal interact with each other 
via a pair-potential we calculate barriers for hopping and exchange diffusion moves over 
island step edges. In this investigation the dependence of the barriers on the island 
size, the misfit and the potential-type is of special interest. We show that the exchange 
barriers are particularly sensitive to these parameters. For example, in the case of large 
misfits the exchange diffusion barrier decreases dramatically and the exchange process 
becomes the dominant diffusion mechanism at island edges. 

In chapter 3, we introduce an algorithm for the simulation of heteroepitaxial growth. 
In order to account for strain effects caused by the atomic mismatch in the crystal it 
is essential to allow for continuous particle positions in our simulations. The main idea 
of our off-lattice method is to compute the barriers for each diffusion event from a 
pair-potential and to use the barriers in a rejection-free Kinetic Monte Carlo (KMC) 
simulation. The main advantage of the KMC method is the feasibility of the relevant 
time scales necessary for the simulation of crystal growth, as realized, e.g., in molecular 
beam epitaxy (MBE). We discuss in this chapter the calculation of activation barriers, 
the relaxation of the crystal and an efficient implementation of the algorithm. 

In the following chapters we discuss the application of our off-lattice KMC method 
to the simulation of three important strain relaxation mechanisms. In chapter 4 we 
investigate the strain relaxation by introduction of misfit dislocations in the adsorbate 
film. In the first part of the chapter we discuss in detail different formation mechanisms 
of dislocations and investigate the so-called critical thickness for the first appearance 
of misfit dislocations. In agreement with various experimental studies the dependency 
of the critical thickness on the misfit is given by a power-law. In the second part 
of the chapter we treat the pseudomorphic region of heteroepitaxial growth and the 
subsequent gradual relaxation of the adsorbate film. These studies are motivated by 
a new in-situ method for the determination of the vertical lattice constant during 
MBE growth. Our simulation results show an excellent qualitative agreement with 
experimental results for ZnSe/GaAs heteroepitaxy. We are able to demonstrate that 
the region of pseudomorphic growth scales with the misfit. 

In chapter 5 follows the investigation of another prominent relaxation mechanism 
in heteroepitaxial growth: the so-called Stranski-Krastanov growth mode, where 3d 
islands self-assemble on a thin adsorbate wetting-layer. We are able to show, that 
an increased barrier for diffusion of adsorbate particles on the substrate is a possible 
kinetic reason for the formation of a stable wetting-layer. The formation of adsorbate 
islands is due to a partial relaxation of adsorbate material, which causes two kinetic 



effects: diffusion on the relaxed island surface is slower than on the pseudomorphic 
strained wetting-layer and a diffusion bias drives particles on the island surface to the 
island center. We further measure the dependence of island size and island density on 
the misfit, the temperature and the particle flux and find good qualitative agreement 
with metal-organic vapor phase epitaxy (MOVPE) experiments. 

Finally, in chapter 6 we consider surface alloying as a possible strain relaxation 
mechanism in multi-component systems. By means of equilibrium simulations we are 
able to show that the competition between binding and strain energy yields the forma- 
tion of regular stripe patterns on the crystal surface. The pattern formation is then 
investigated in a kinetic model: depending on temperature, misfit and potential-type 
we find as well the pattern formation as the experimental reported ramification of the 
structures. The comparison with a lattice gas model shows that the observed pat- 
tern formation can solely be due to kinetic effects but the misfit is essential for the 
ramification and stabilization of the surface structures. 

In the appendix of the work we discuss origin and properties of the used pair- 
potentials. 
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Chapter 1 



Theoretical descriptions and models 
of heteroepitaxial growth 

In recent years heteroepitaxial growth has been a field of intense study. This is mainly 
due to the fact that countless technical applications - including laser diodes, solar cells 
and magnetic or optic data storage devices - are heterosystems. In this work we focus 
in particular on atomic size mismatched systems where the participating materials 
crystallize in the same lattice structure but have different lattice constants. 

In technical application this so-called misfit can lead to both wanted and unwanted 
effects. On the one hand one often wishes to deposit a smooth adsorbate layer of certain 
thickness on a given substrate with a different lattice constant. This is for example 
needed in the fabrication of computer memory chips. Here, the relaxation of strain 
due to the misfit can cause perturbations of the lattice structure (e.g. dislocations) 
or modulations of the surface which may affect the characteristics of the device in a 
negative way. 

On the other hand a moderate misfit can lead to self-assembly of islands which is an 
interesting process for the fabrication of so-called quantum dots. These quantum dots 
are solid state structures typically made of metals or semiconductors which confine a 
small number of electrons to a small space used for optoelectronic devices or single 
electron transistors. However, in both cases the understanding and control of misfit 
caused phenomena is essential. 

1 .1 Microscopic processes on crystal surfaces 

A particularly important technique for the fabrication of heterostructures is the molec- 
ular beam epitaxy (MBE) and related techniques: in an ultra high vacuum (UHV) 
environment the substrate surface is exposed to a uniform flux F of adsorbate parti- 
cles, which are evaporated in a thermal effusion cell. The interplay of three relevant 
microscopic processes on the crystal surface determines the morphology of the growing 
adsorbate film: deposition, desorption and surface diffusion (see fig. 1.1) of adsorbate 
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particles [1-3]. 
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Figure 1.1: Schema- 
tic depiction of rele- 
vant processes during 
MBE growth. 



1.1.1 Deposition 



From the vapor beam adsorbate particles arrive with a flux F - normally measured in 
monolayers per second (ML/s) - at a random position on the surface of the crystal. 
Since the temperature of the effusion cell is typically higher than the substrate tem- 
perature, the energy of the newly deposited particle is higher than that of particles 
already in thermal contact with the surface. This may result in a increased mobility of 
the arriving particles, which is accommodated in some simulation models by effective 
rules like incorporation or downhill funneling (see e.g. [4]). In the latter case the initial 
motion of the particle is biased to energetically favorable sites with a high coordina- 
tion number where the particle then stick to the surface. After the deposition process 
has ended the particle is considered to be in a chemisorbed state at the bottom of a 
potential well - the so-called binding state - thermalized at the crystal temperature. 

In the following we consider only the chemisorbed state and neglect physisorbed 
states, in which the adsorbed particle is held at the surface by much weaker Van der 
Waals forces. This can lead to an enhanced mobility of the adatoms and is addressed 
in detail in recent publications [4,5]. 

1.1.2 Desorption 

A competing effect of deposition is the desorption of a bound adsorbate particle. The 
probability of desorption depends on the depth of the potential well - the so-called 
binding energy Eb - and the temperature T of the crystal surface. Thermal fluctuations 
tend to drive the adsorbate particle back to the gas phase with a rate proportional to 
exp(—Eb/kT). The desorption rate Rj is given according to an Arrhenius law [2]: 



Here Uq is the attempt frequency which is on the order of magnitude of the Debye 
frequency of the crystal. The binding energy E b depends both on the specific particle 
types and the local geometry of the surface i.e. the coordination number of the binding 
site. In heteroepitaxial growth E b depends also on the misfit between adsorbate and 
substrate particles and varies locally depending on the local strain. 



Rd = ^oe kT ■ 



(1.1) 
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Due to the high binding energies (typically E b 2.5eV) the desorption rate is small 
compared to the rates of deposition and diffusion. Thus for many material systems 
under typical MBE conditions the desorption of chemisorbed adsorbate particles is 
negligible. 

1.1.3 Surface diffusion 

The surface diffusion of adsorbate particles is described by an activated process where 
particles jump laterally along the surface from one binding place to another. In this 
process they have to overcome an energy barrier, the so-called activation energy E a . 
At the transition state - the transition energy E t - the particle is less well bound than 
at the binding site but still far from zero of energy (corresponding to desorption). The 
activation energy is therefore given by 

E a = E t -E b (1.2) 

and again according to an Arrhenius law the rate for surface diffusion R results to 

R = UQ e~Tf. (1.3) 

In the multi-dimensional case (d > 2) the transition state corresponds to a first order 
saddle point in the potential energy surface (PES) [6,7] with a maximum in the direction 
of diffusion and minimum along all other coordinates of the crystal surface [8]. The 
binding site is represented by a minimum in the energy landscape (see fig. 1.2). 




Figure 1.2: PES for a test par- 
ticle on a plain surface. The par- 
ticles of the crystal interact via a 
3d cubic Lennard- Jones potential 
(cf. A). 



Potential energy surface 

Consider a crystal containing of n — 1 particles interacting with each other via a pair- 
potential. In order to compute the PES of the crystal surface a test particle n is moved 
in small steps across the surface. After each step the total potential energy of the 
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n-particle system is minimized by variation of all particle coordinates plus the test 
particle's coordinate perpendicular to the surface. Potential energy and coordinates 
of particle n result then in the PES. This so-called Molecular Static method for the 
calculation of the energy surface is explained in detail in chapter 2. The energy surface 
represents the stable (minima) and metastable (saddle points) sites for an adatom at 
crystal temperature T = (disregarding the vibrations of the lattice). Figure 1.2 shows 
such a PES for a particle on a plain 2d surface. 



Diffusion bias 

Figure 1.3 shows the PES of a 2d Lennard-Jones crystal, where due to the isotropy 
of the pair-potential the particles arrange into a triangular lattice. On the substrate 
with lattice constant a s a monolayer island of adsorbate particles (lattice constant a a ) 
is placed. The test particle moves from the middle of the island towards the right side 
of the crystal. As one would expect for a triangular lattice the minima of the PES 
correspond to the sites in the middle between two particles in the underlying layer 
(bridge sites). As a 2d crystal is considered here the transition states are represented 
by maxima in the PES and coincide with the vertices of underlying particles (top sites) . 
When the test particle crosses the edge of the island it has to overcome an increased 
energy barrier E s . This so-called Schwoebel barrier [9, 10] results from the reduced 
number of binding partners at the rim of the island. The probability of being reflected 
at the edges is therefore higher than the one for jumping from the island. On the 
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Figure 1.3: PES for a test 
particle on the surface of a 2d 
Lennard-Jones crystal. The 
misfit is e = 5%. Note that 
the the diffusion barriers rise 
from the island center to the 
edge. 



other hand a test particle moving towards the bottom of a monolayer islands has to 
overcome a decreased barrier in order to reach the island step. Within the range of 
the interaction potential between the particles of the system islands act attractive to 
surrounding particles. Because of the increased number of binding partners a particle 
located at the step edge is much stronger bound than on a flat surface. 



S 



Another phenomenon of biased diffusion on top of islands is especially pronounced 
in heteroepitaxial growth, when the misfit 



a s 

between substrate and adsorbate becomes nonzero. In order to minimize the potential 
energy of the system, adsorbate islands always try to achieve their favored lattice con- 
stant a a . Because of the higher mobility this works most efficiently for particles at the 
edges of islands, whereas the center of an island is still strained and the lattice constant 
there is closer to a s . As a consequence of the inhomogeneous relaxation the binding 
and transition energies for diffusion depend both on the position of the adatom on the 
island and the misfit between adsorbate and substrate. The diffusion close to the edge 
is different from that near the island center [11]. For negative misfit (e < 0) this leads 
to a diffusion current from the island center towards the edges. This can e.g. result in 
a reduction of the next layer nucleation rate. On the other hand in the case of positive 
misfit (e > 0) the diffusion is biased to the island center and next layer nucleation can 
be enhanced. Biased diffusion from the island edge to the center can also be concluded 
from figure 1.3, where the misfit between adsorbate and substrate is e = 5%. 

In systems with pair-potentials such a behavior of the diffusion barriers is a general 
phenomenon [11] and is observed for all types of pair-potentials used in this work: if 
the diffusing particle feels a smaller underlying lattice constant than its natural one 
(i.e. diffusion on a compressive strained crystal) the diffusion barrier is decreased and 
diffusion becomes faster. For diffusion on a tensile strained crystal the situation is con- 
trariwise and the activation energy is increased. This can be understood in an intuitive 
way: compression of the lattice moves the diffusing particle a bit away from the surface. 
For that reason the particle experiences a less undulated PES. For the extreme situation 
of a very large compressive strain the underlying crystal can be viewed as a continuous 
film with no discrete binding sites and the diffusion barriers therefore vanish. In the 
limit of large tensile strain diffusion is equivalent to breaking a pair of atoms and build 
a new one resulting in a diffusion barrier equal to the pair binding energy. As shown 
in [11] for the Lennard- Jones potential the dependence between activation energy and 
strain is linear, at least for moderate values of e. Given real materials metallic systems 
show the same trend of the strain dependence. In case of semiconductors the strain 
dependence of the diffusion barriers cannot be explained that easily [12]. For exam- 
ple first -principles calculations for the In/GaAs(001) surface showed that the diffusion 
barrier has a non-monotonic strain dependence with a maximum at compressive strain 
values [13]. 

1.1.4 Exchange diffusion 

Although the above described hopping diffusion is the intuitive way one would think of 
atoms moving over an surface, there is another important diffusion mechanism. Espe- 
cially in the case of face-centered-cubic (fee) metals the so-called exchange diffusion 
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becomes relevant. In the latter, the adatom takes the place of a lattice atom whereas 
the displaced atom becomes the new adatom and continues the diffusion [14]. This 
mechanism is most likely for diffusion over descending steps where the smaller number 
of neighboring atoms makes the displacement of an atom from the step rather easy (see 
chapter 2). 

1.2 Strain relaxation mechanisms 

In the following the bulk lattice constant of the substrate and the adsorbate will be 
denoted as a s and a a . The misfit between the substrate and the adsorbate film is given 
according to equation (1.4). 

If the misfit is not too high << 1) the adsorbate is coherent with the substrate 
during the early stages of growth. In this state the crystal topology is that of a perfect 
crystal, i. e. each particle has the same coordination number and its nearest and next- 
nearest neighbors form the same geometrical figure with only slightly modified distances 
[3,15]. 

As the thickness h of the adsorbate film increases the elastic energy of the film rises 
and there are two possible relaxation mechanisms: the introduction of misfit dislocations 
in the adsorbate layer and the formation of three-dimensional adsorbate island. 

1.2.1 Misfit dislocations 

An important strain relaxation mechanism which is often ineluctable for large adsorbate 
layers in large misfit heteroepitaxy is the formation of dislocations, where the strain 
energy is released by plastic deformation. In this case, if the absolute value of the 
misfit between substrate and adsorbate |e| is sufficiently small [15], the crystal topology 
is only weakly perturbed in large domains separated by lines. Along these lines, called 
misfit dislocations the perturbation is large. The adsorbate film thickness at which 
dislocations first occur is known as critical thickness h d c and qualitatively increases with 
decreasing values of 

The existence of a dislocation is indicated by a non-vanishing dislocation-displace- 
ment vector, the so-called Burgers vector b [16]. As a result of continuous elasticity 
theory [15,17,18] only the component of b parallel to the substrate/adsorbate interface 
contributes to the relaxation of strain. With the help of the Burgers vector dislocations 
can be categorized in climb and glide dislocations. Climb dislocations (fig. 1.4(a)) - 
with a Burgers vector parallel to the interface - relax the elastic energy best but the 
system has to overcome a high activation energy to create them (see e.g. [19]). The 
easier to form glide dislocations (fig. 1.4(b)) have a component of b vertical to the 
interface, which does not contribute to the relaxation. 

If the crystal topology is only perturbed near the dislocation line and far from the 
interface the topology of the crystal is the same as in the coherent state the Burgers 
vector is a lattice vector and the dislocation is called perfect (fig. 1.4(a),(b)). Otherwise 
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(a) e = 10% (b) e = 5% (c) e = 6% 

Figure 1.4: Sections of crystals obtained in simulations with perfect climb (a) and 
glide (b) dislocations. The dashed arrows show the Burgers circuit, which is drawn to 
determine the dislocation Burgers vector b. Panel (c) shows a partial glide dislocation 
with an partial atomic step at the surface. 

if h is not large enough partial dislocations are formed (fig. 1.4(c)), characterized by 
a Burgers vector which is a rational fraction of the lattice vector. In this case lattice 
planes have a discontinuity when crossing the glide plane and a partial atomic step 
appears on the crystal surface. 

1.2.2 Island formation 

Mismatched epitaxial films relax their strain not only by the introduction of misfit 
dislocations, but also in an elastic way by deformation of the surface or the so-called 
2d — 3d transition. Historically, three growth modes are distinguished in heteroepitaxial 
growth: the Frank-Van der Merwe (FM), the Volmer-Weber (VW) and the Stranski- 
Krastanov (SK) type of growth. Which of these growth modes in thermal equilibrium 
occurs depends upon the relative magnitudes of the surface energies j s , 7 a and 7; (h) of 
the substrate, the adsorbate film and the interfacial energy, respectively. 7 S and 7 a are 
the values for the semi-infinite crystals. The strain energy depending on the thickness 
of the adsorbate film h has been absorbed in 7* here [20,21]. 

Frank-Van der Merwe growth mode 

In the FM growth mode the adsorbate forms a flat film on the substrate in a layer-by- 
layer way. It occurs when 

A 7 = 7 a + 7i(/i)-7 s <0 (1.5) 

for all film thicknesses h. It is observed for several material systems with small mis- 
fits |e| < 2% at low deposition fluxes and high temperatures, where 3d nucleation is 
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suppressed and adsorbate particles are mobile enough to reach surface steps. Since the 
strain in the adsorbate film rises with increasing thickness h the FM growth mode is 
metastable [22]. It can be relaxed mainly by the introduction of misfit dislocations or 
the deformation of the surface, known as Asaro-Tiller-Grinfeld instability [15]. 

In the latter case a weak perturbation of the adsorbate surface leads to a quasi- 
periodic modulation: For instance, in a heteroepitaxial system with a small positive 
misfit the stress is relaxed at the peaks of a weak surface modulation, because the film 
is tensed there and the adsorbate can get closer to its own lattice constant a a . On 
the other hand stress is concentrated at the valley, where the adsorbate film is under 
compression. This difference in strain energy density causes mass transport by surface 

Figure 1.5: Strain-induced epilayer 
roughening: strain is reduced at the peak 
of the modulation and increased at the 
troughs. The arrows symbolize the direc- 
tion of surface diffusion. 

diffusion from high to low strained regions (see fig. 1.5) and leads to a quasi-periodic 
three-dimensional modulation of the surface of a rather thick adsorbate film [23]. For 
instance, this thickness is several hundreds of monolayers for Sio.84Geo.i6 o n Si(OOl) with 
a misfit of e « 0.6% [24]. 

Stranski-Krastanov growth mode 

If the FM condition (1.5) is only satisfied for a small number of adsorbate layers SK 
growth is energetically possible. In this growth mode coherent three-dimensional (3d) 
islands form on a so-called wetting-layer of coherently growing adsorbate (see fig. 1.6). 
Typically the thickness of this wetting-layer is between one and four monolayers. The 
formation of SK islands can be understood as a phase transition: the 2c? wetting- 
layer is metastable and grows up to a supercritical thickness h*, when more stable 
islands begin to form. Further growth of these islands is fed both by capturing of newly 
deposited adsorbate particles and the decomposition of the supercritical layer. After the 
transition, the thickness of the wetting-layer decreases to a stationary value, commonly 
denoted by h c . It is important to stress, that this process even takes place without any 
further deposition of particles [25]. 

Figure 1.6: Schematic representation 
of the Stranski-Krastanov growth mode, 
where 3d adsorbate islands form on an ad- 
sorbate wetting-layer. 

The exact mechanism of island formation is still under discussion. It seems clear, 
that monolayer (2d) islands - located on the wetting-layer - play an important role as 
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precursors for the formation of 3c? islands. One theoretical model is that at a critical size 
these 2d islands become unstable. Because of the strain the more weakly bound edge 
atoms jump to the next island level and build up the second island layer. This process 
is repeated in the transformation from bi- to three-layer islands, and so on [26,27]. 
It is known as the spontaneous 2d — 3d transition which results in the formation of 
self-assembled quantum dots. 

This transition involves an activation barrier E2d-3d which the system has to over- 
come. Otherwise the film thickness may increase until the introduction of dislocations 
becomes favorable. Figure 1.7 gives schematically the energetic situation for the two 
cases of SK and dislocation-relaxed FM growth. First the total energy of the system 



introduction of misfit 
dislocations 




Figure 1.7: Total energy of a 
strained heteroepitaxial system ver- 
sus adsorbate layer thickness h. 3d 
island nucleation starts at h* c) dislo- 
cations nucleate at h d c [22]. 



nucleation of 3d islands 



3 

h/ML 



decreases due to the energy contribution from the substrate/adsorbate interface until 
the substrate is covered by a complete monolayer of adsorbate. Then the elastic energy 
in the strained adsorbate film increases linearly with each film layer. If the systems 
does not manage to overcome the barrier E 2 d~3d FM growth is maintained until the 
introduction of misfit dislocations at a film thickness h d c . If on the other hand growth 
conditions favor the SK mode 3d island nucleation starts at h*. This SK growth is 
observed in various strained heteroepitaxial systems of the IV-IV, III-V and II- VI 
families of semiconductors (for a review see [22]). In all these cases the misfit is positive 
and quite large (2% < e < 7%). 



Volmer-Weber growth mode 

In the case of VM growth equation (1.5) is not fulfilled and 3d adsorbate islands form 
directly on the substrate (see fig. 1.8). Due to the fact that there is no technological 
application of incoherent islands - where misfit dislocations are introduced to relieve 
the strain energy - coherent VM growth is of special interest. In this case the base of 
the islands is still constrained by the substrate but the adsorbate can reach its own 
lattice constant a a on the top, the sides and to some degree in the island center. The 
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mechanism of island formation should be similar to the process of 2d — 3d transition 
described above. 




Figure 1.8: Schematic representation of 
the Volmer- Weber growth mode, where 3d 
adsorbate islands form directly on the sub- 
strate. 



Coherent VW growth is observed mainly in heteroepitaxial systems with large pos- 
itive misfit, e.g. for ZnTe/ZnSe [28], Mn/Si(lll) [29] or Si/Ge(lll) [30]. 

1.2.3 Surface confined alloying 

A further strain relaxation mechanism, which arises generically in systems dominated by 
atomic size mismatch is surface confined alloying of the materials [31]. Surface alloying 
is observed for both cases: two component systems with mixing of adsorbate and sub- 
strate particles (e.g. Na/Al(lll), K/A1(111) [32,33], Au/Ni(110) [34], Ag/Pt(lll) [35], 
Sb/Ag(lll) [36]) and three component growth with alloying of two adsorbate species on 
a different substrate (e.g. CoAg/Ru(0001) [37-40], CoAg/Mo(110) [41], FeAg/Mo(110) 
[41], AgCu/Ru(0001) [42], PdAu/Ru(0001) [43]). 

In principle, e.g. for the case of a binary system, surface alloying is again understood 
in terms of surface and interface energies [35]. If the interfacial energy 7, < adsorbate 
and substrate can lower their energy by intermixing. Otherwise the adsorbate material 
will segregate. 



Figure 1.9: Schematic representation of 
surface confined alloying of two adsorbate 
species (dark and light gray) on a different 
substrate. 

Also for the three component system the strain relaxation due to surface alloying 
can be understood intuitively. Consider the situation of figure 1.9, with two adsorbate 
species A, B which imply a misfit of the same absolute value but opposite sign with 
the substrate. As long as no further difference between the adsorbate particles exists 
and in particular the binding E AB between the two species is the same as for the A- A 
and the B B interaction, an alternating arrangement of both adsorbate types is likely 
to be the energetically most favorable state. A weaker A-B interaction can complicate 
the situation and causes a competition between alloying and the introduction of misfit 
dislocations as the preferred relaxation mechanism. 
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1 .3 Models of heteroepitaxial growth 



Various methods have been proposed recently to model strain relaxation and growth of 
atomic mismatched systems by means of computer calculations and simulations. Due to 
the steady growth of computer power they became an important tool for analyzing and 
understanding microscopic processes and their effects. In this section we give a short 
overview of methods for the computational treatment of, in particular, heteroepitaxial 
growth. 

1.3.1 Density-functional-theory 

For calculations and simulations of growth processes it is desirable to describe elemen- 
tary processes, like the interaction between atoms or molecules, on a high quality level 
i.e. without making empirical assumptions. In contrast to classical approaches de- 
scribing the atom-atom interaction using empirical potentials or simple bond-counting 
rules, density-functional-theory (DFT) (see e.g. [7]) accounts for the quantum mechan- 
ical nature of the electrons. 

For DFT calculations only very few assumptions about the electronic structure of 
a poly-atomic system have to be made: an approximate density functional is used 
to solve separate, Schrodinger-like equations for all electrons. Electron many particle 
(exchange or correlation) effects are added in terms of an additional potential. 

With respect to the high computational effort of the method it is not possible to 
use DFT for in-situ calculation of e.g. diffusion barriers in Monte Carlo simulations. 
Only the calculation of single events like e.g. the adsorption of an As 2 molecule on a 
GaAs surface at zero temperature and pressure is within the scope of the method right 
now [7]. 

However, DFT methods were used with some success for energy calculations of ex- 
emplary situations like found in heteroepitaxial growth. For example in [44, 45] the 
total energy of a system consisting of a pseudomorphic InAs layer located on a GaAs 
substrate was compared to the energy of a system with an InAs island of certain shape 
and size placed on a thinner InAs film. DFT methods were used for the calculation 
of surface energies, whereas the elastic energy was treated by means of classical elas- 
tic continuum equations. The investigation yielded results about energetically most 
favorable island shapes and sizes and predicted a non-vanishing wetting-layer due to 
energetic reasons for the InAs/GaAs system. 

In conclusion, though DFT is at present not suitable for large scale simulations 
of heteroepitaxial growth processes, it is an appropriate method for material specific 
calculations regarding prototype situation in atomic mismatch systems. 

1.3.2 Molecular Dynamics simulations 

Since the more or less exact calculations of chemical bonds between atoms or molecules 
are computationally far too demanding for the simulation of many particle systems 
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and relevant time scales often empirical potentials are chosen. These potentials can 
give the interaction strength between particles of the systems as a function of e.g. 
the particle distance and the bond directions. Material specific potentials are fitted 
in order to reproduce bulk properties like elastic constants, the vacancy formation 
energy, the stacking fault energy, surface energy, and phonon frequencies [46]. One 
prominent example for empirical material specific potentials are the Embedded Atom 
Method (EAM) potentials used for the modeling of metallic systems. EAM potentials 
are calculated as a sum of pairwise interactions between the particles of the system and 
a many body term. 

Empirical potentials for the simulation of semiconductor properties are more com- 
plicated and computationally more demanding. An important set of potentials for semi- 
conductor systems are e.g. the so-called Tersoff-potentials [46,47]. Tersoff-potentials 
are based on the concept of bond order: the strength of a bond between two atoms is 
not constant, but depends on the local environment. The basic idea is that the bond 
between two atoms i and j is weakened by the presence of other bonds i — k also involv- 
ing the atom i. The amount of weakening is determined by where these other bonds 
are placed. Hence angular terms become necessary for the construction of a realistic 
model. 

Given at least an empirical potential for the interaction between particles of the 
system, Molecular Dynamics (MD) simulations are clearly the most realistic and de- 
sirable method for the simulation of heteroepitaxial growth. In MD simulations the 
time evolution of a particle system is described by integrating the equations of motion 
using the interaction potential for the computation of the forces on each particle. MD 
techniques are applied with success to e.g. the investigation of dislocation formation 
and motion [48-50] or the formation of grain boundaries [51] for rather large system 
sizes (up to 10 6 particles). 

However, MD simulation suffer generally from the restriction to short physical times 
on the order of 10~ 6 s or less. The simulation of crystal growth - like e.g. in the MBE 
environment - requires the coverage of seconds up to minutes. That is because the 
properties of a growing crystal are ruled by rather rare thermally activated processes 
like surface diffusion jumps from one local minimum to a neighboring one. These events 
normally occur with an exponentially decreasing probability (see eq. (1.3)) [7]. 

MD simulations are altogether surely an adequate method for the investigation of 
concerted moves involving many particles like dislocation motions or interdiffusion pro- 
cesses at the interface between two material types. But MD methods are not suitable 
for the simulation of growth processes since relevant time scales are currently not fea- 
sible. The simulation of heteroepitaxial growth requires a further simplification with 
concentration on the relevant rare events. 

1.3.3 Off-lattice Monte Carlo simulations 

Monte Carlo simulations have been proved to be a powerful tool for the simulation of ho- 
moepitaxial growth (see e.g. [4]). Here rare events (e.g. diffusion events) are performed 
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according to their probability, given by an Arrhenius law (see eq. (1.3)). However, for 
the simulation of heteroepitaxial growth it is essential to allow for continuous particle 
distances in order to account for strain effects in the crystal. 

The computationally least demanding simulation technique which fulfills this re- 
quirement is the so-called ball and spring model [52-57] (see fig. 1.10). Here the ac- 




Figure 1.10: Schematic representation of 
the ball and spring model. The elastic en- 
ergy is determined by the spring constant 
k x and the natural length a x . 



tivation energy for a diffusion jump of a surface particle is split into bond and strain 
energy E a = E^ond — E stra in- The bond energy is determined by the exact number of 
nearest (nn) and next-nearest neighbors (nnn) e.g. by a simple bond-counting rule. 
The strain energy for a site i is obtained by the difference of the system's elastic en- 
ergies with site i occupied and unoccupied. The elastic energy is given by harmonic 
interactions (springs) between an atom and its nn and nnn. Interaction between sub- 
strate particles is represented by the spring constant k s and the natural length a s . The 
adsorbate-adsorbate interaction is given by k a , a a — (1 + e)a s , accordingly. 

Since the calculation of E bond and E strain requires the number and positions of sur- 
rounding particles (nn and nnn) of each particle this method does not allow for the 
simulation of misfit dislocations. The method also suffers form the rather rough de- 
scription of the elastic energy by simple harmonic interactions. However, ball and spring 
models yield some interesting results on the growth of coherent islands and the different 
growth modes. 

A step towards a more realistic treatment of the binding between particles are 
continuous-space Monte Carlo simulations [58-60]. The interaction between the par- 
ticles of the system is here given by pair-potentials. A surface move of a particle is 
e.g. accepted with the probability exp (— E^/kT), where Ef, is the binding energy of 
the particle due to the interaction with surrounding particles of the crystal within a 
certain range. Since - unlike in the ball and spring simulations - number and positions 
of surrounding particles are not required for the calculation of Eb this method allows 
for the simulation of dislocations [60]. 

In this work we go a different way and put forward an off-lattice Kinetic Monte 
Carlo (KMC) method introduced in [61]. In the KMC method the activation energy 
for a hopping diffusion jump from one local minimum to another is calculated with 
respect to a pair-potential that mediates the interaction between the particles of the 
system. The method will be explained in detail in chapter 3. In chapter 4 we analyze 
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the formation of misfit dislocations with this method. Chapter 5 gives investigations 
of the SK growth mode and in chapter 6 we use the off-lattice KMC method for the 
simulation of multicomponent growth. 
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Chapter 2 

Diffusion barriers on step edges 



In this chapter we focus on the influence of the misfit between substrate and adsorbate 
on the diffusion barriers across steps. 

The barrier for the downward movement from the top of an island plays an important 
role for the growth mode: high barriers for downward movement can lead to rough 
growth, whereas a low barrier for descending jumps favors smoother layers at a given 
temperature. 

Large effort has been spent on the calculation of diffusion barriers for various mate- 
rial systems by the means of empirical potential methods [62-66] or density functional 
theory (DFT) [67,68], mostly for the case of self-diffusion. But to our knowledge there 
is no systematical study on the influence of the lattice misfit between substrate and the 
adsorbate island on the downward barriers. Since the control of the surface morphology 
is a major goal in MBE techniques the influence of the misfit on barriers for descending 
diffusion moves is of special interest. 

By means of a Molecular Static method we study these barriers for different misfits, 
island sizes and interaction potentials for the 2d and 3d case. The difference of the 
barriers for hopping and exchange diffusion are of particular interest. The computation 
of these barriers is also of some relevance for the following parts of this work: in the 
case of off-lattice simulations it is important to know whether the rather complicated 
concerted moves - like e.g. exchange diffusion - should be taken into account. 

2.1 Hopping and exchange diffusion 

As mentioned in chapter 1 surface diffusion of an adatom can proceed by two different 
processes: 

• Hopping diffusion: The adatom jumps from one minimum in the potential energy 
surface (PES) to another by overcoming the activation barrier E a = E t — E b , 
where Ej, is the energy of the binding state and E t is the energy at the transition 
state. 
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• Exchange diffusion: The adatom replaces a surface atom and this surface atom 
resumes further diffusion [14]. The activation energy of the exchange process can 
also be written as E a = E t — E\,. 

The latter diffusion mechanism may become particularly important at descending steps. 
Due to their lower coordination number atoms at the step edge can be more easily 
pushed away by the adatom (see fig. 2.1), whereas for the case of hopping over the rim 

Figure 2.1: Schematic representation of a ex- 
change diffusion move at a descending step edge 
for the 2d case. The adatom takes the place of the 
surface particle at the island's edge. 

of an island in many material systems the adatom has to overcome an energetically 
unfavorable situation due to the weaker binding during the diffusion step. This is 
especially pronounced in the 2d case (see fig. 2.2). 

Figure 2.2: Schematic representation of a hop- 
ping diffusion move over the rim of an island. 



2.2 Calculation of diffusion barriers 

The technique we apply here in order to compute the barriers for downward hopping and 
exchange diffusion at a descending step edge is the so-called Molecular Static method 
(see e.g. [61-66]). 

Consider a crystal with n—1 particles interacting via a pair-potential Uij which is a 
function of the distance \fij\ between two particles % and j of the crystal. An additional 
test particle n is placed on the crystal's surface and the system's total potential energy 

n— 1 n 

E tot = E E u v 

8=1 j=i+l 

is minimized by the means of a conjugate gradient method [69]. The coordinates 
of the n — 1 crystal particles and the test particle's coordinate perpendicular to the 
surface are varied. After reaching the minimum energy configuration of the system E tot 
together with the position of the test particle is noted and the test particle is moved 
by a small step 5 in a given direction. This procedure is repeated until the particle 
n has reached a pre-determined position. The recorded particle coordinates together 
with E tot then result in the PES of the crystal (see figures 1.3 and 1.2 for the 2d and 
3d case, respectively). 
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2.3 Calculations for the two-dimensional case 



For the calculation of the downward hopping barrier in the 2d case the test particle is 
set on top of the island, near the step edge (see fig. 2.2) and moved towards the step 
edge. The activation energy E ajhap is given according to equation (1.2) (see fig. 2.3). 
For the exchange diffusion move the test particle is drawn out of the step edge and the 
particle on top of the island relaxes into the resulting gap (see fig. 2.1). Figure 2.3 
shows the resulting exchange barrier E aex . 



Figure 2.3: PES for hopping and 
exchange diffusion from the top of 
an island in two dimensions. The 
particles of the system interact via 
a Morse potential (a = 5.0 and 
e = 0). Note that in case of ex- 
change diffusion x gives the posi- 
tion of the edge particle, whereas 
in case of hopping diffusion x de- 
notes the position of the hopping 
particle. 




2.3.1 Parameters 



For the calculation of the PES we take a 15 layers thick substrate (h=15), which contains 
L = 70 particles in each layer. Periodic boundary conditions are applied in the lateral 
x-direction. To stabilize the crystal during the relaxations, the bottom layer is frozen, 
that means particles in this layer are not allowed to relax. 

On the substrate an adsorbate monolayer island is placed consisting of 3 < I < 20 
particles. As interaction potential for the particles of the crystal the Lennard- Jones 
12, 6 potential 
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(2.2) 



(2.3) 



and Morse potentials 

for values of a between a = 4.5 and a = 7.0 are chosen (also see appendix A). 

In the following the potential depth is set to Eij = l.OeV for all types of interactions. 
The potential is cut off at a distance r cut = 6.0r , which is perfectly justified since all 
used pair-potentials converge fast to = for the particle distance — > +oo. The 
results were also compared to preliminary calculations with r cut = 4.0r and r cut = 
12.0ro which yield no significant differences in the obtained barriers. 
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The misfit between the substrate and the island is varied between e = — 14% and 
e = 14% in steps of 1%. For each value of the island size / and the misfit e the PES 
near the step edge is calculated for both possible diffusion mechanisms. The barriers 
E ay hop and E a ^ ex then result from the obtained energy landscapes. 



2.3.2 Exchange vs. hopping 

In the following we examine the difference between the barriers for exchange and hop- 
ping E a y lop — E aex as functions of island size and misfit. Figures 2.4(a) and 2.4(b) show 
this difference for a Morse potential (a = 6.0) as a function of I. Each curve represents 
a value of e. The characteristics of these curves are typical for all used potentials. First 



> 

CD 



LU 

I 



LU 



0.2 r- 

0.1 
0- 
-0.1 
-0.2 
-0.3 
-0.4- 
-0.5 
-0.6- 
-0.7- 



(a) 



- -c © - - 

- -O G- — 



10 



12 



> 

CD 



LU 

I 



LU 



1 

0.8 
0.6 
0.4 
0.2 


-0.2 
-0.4 
-0.6 
-0.8 
-1 



(b) 



o- -o- -o — o- - 



^_ -o- -o- -» -°- 

_o_ -o- -o- -o — 0- -o- 
1 -e- -o- -o- -o- -o- -o- 



10 12 
I 



14 



16 



18 20 



Figure 2.4: E a ^ hop — E afiX as a function of island size / for different misfits e. (a) from 
down to top: e = —0.01 to e — —0.14 in steps of 0.01, (b) from down to top: e = 0.0 
to e = 0.14 in steps of 0.01. 



E a) hop — E a ^ ex increases with / until the difference becomes more or less constant. If \e\ 
exceeds a critical value e c the difference becomes positive, which means E a ^ ex < E a j lop 
and thus exchange diffusion becomes likely. The absolute value of the critical misfit 
varies for positive and negative misfit. For each \e\ > \e c \ there is a critical island size 
l c : for island sizes / > l c the barrier for exchange is smaller than the barrier for hopping. 

This behavior is easy to understand if one considers the position of the particle which 
is drawn out of the island edge. For the homoepitaxy situation (e = 0) all particles of 
an island - and therefore also the edge particles - match perfectly the lattice structure of 
the substrate. But at a given island size with increasing misfit the edge particles move 
toward the top site of an underlying substrate particle and become less well bound. 
The same is true for a given misfit and increasing island size: the larger the island the 
less favorable is the position of the edge particles. 

This affects the exchange diffusion more than the hopping diffusion: in the hopping 
diffusion mode the active particle always has to overcome an energetically unfavorable 
top position, whereas in exchange diffusion a weakly bound edge particle is much easier 
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kicked out of its place. In case of e.g. the a = 6.0 Morse potential shown in figure 2.4 
at e = 0.12 and I = 9 the hopping barrier has 84% of its value at I = 3, whereas the 
exchange barrier dropped to 28% of its value at I — 3. 

To explain the mechanism of exchange diffusion in more detail figure 2.5 shows a 
PES for a high positive misfit (e = 0.12) at island size I = 8 in case of the Lennard- 
Jones 12, 6 interaction. In figure 2.6 corresponding snapshots of the particle positions 
during the exchange process are shown: 




Figure 2.6: The particle positions during the exchange process corresponding to the 
PES given in figure 2.5. Further explanations are given in the text. 



(T) Figure 2.6 (T) gives the positions of the particles before the exchange diffusion 
process starts. This state corresponds to a local minimum of the PES. Note 
that the edge particle on the right hand site of the island is nearly on top of an 
underlying substrate particle. 
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(2) The edge particle is drawn over the top of the underlying substrate particle indi- 
cated by a local maximum of E tot . 

(3) Now the edge particle has reached a bridge site between two particles of the 
substrate. Due to the high positive misfit also the particle on top of the island 
runs across a favorable binding position. This leeds to a further local minimum 
in the PES and corresponds to the formation of a dislocation at the edge of the 
monolayer island. 

(4) In the following the edge particle has to overcome another top site resulting in a 
local maximum of E tot . 

(5) Finally the exchange diffusion is finished resulting in an I = 9 monolayer island. 

(B) Due to the high misfit the bond of the edge particle is weak and a further jump 
away from the island's edge involves an only small barrier. 

(7) Because of the spatial proximity of the island the second next binding place away 
from the island results in a rather deep local minimum of the total energy of the 
system. 

We conclude from this investigation on the exchange mechanism that in the large 
misfit and large island regime due to the unfavorable binding position of the edge 
particle exchange diffusion becomes more likely. But one has to consider that at least 
in case of positive misfits the system has to pass through a further binding state, which 
is related to the introduction of a misfit dislocation on the island edge. 

2.3.3 Influence of the potential 

We analyze now the influence of the used potential Uij on the calculated barriers. Since 
the exchange move involves a lot of stretching and compressing of the participating 
particles (see also fig. 2.6), the exchange barrier should be especially susceptible to the 
characteristics of the used interaction around its equilibrium distance. 

This assumption is confirmed by figure 2.7(a). Here the exchange barrier E afiX 
for a island of size / = 7 is shown as a function of the misfit. As one would expect 
exchange diffusion involves the highest barriers close-by e = for all used potentials. 
In comparison to figure 2.7(b) it becomes clear that the steeper the potential around the 
equilibrium distance, the higher becomes the barrier for exchange diffusion. Therefore 
the Morse potential with a = 7.0 gives the highest and the Morse a = 4.5 potential 
the lowest exchange barriers at the same misfit. Of particular interest is the case of 
the Lennard-Jones interaction: as figure 2.7(b) displays the Lennard- Jones potential 
overlaps to a high degree with the a = 5.0 Morse potential for particle distances greater 
than the equilibrium distance. This results in the collapse of the exchange barriers of 
both potentials for large positive values of e. For r\j < r it behaves more like a Morse 
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Figure 2.7: (a) E a ^ ex as a function of the misfit for several Morse and the Lennard- 
Jones 12, 6 potential, (b) The characteristics of the used pair-potential £7^ around their 
equilibrium distance. For a better comparability the minimum of the Lennard- Jones 
potential is shifted to r\j = 1.0. 



a = 6.0 potential which is reflected by the shift to the Morse a = 6.0 results for highly 
negative values of e. 

It is also seen from figure 2.7(a) that for a given absolute value of the misfit \e\ the 
exchange diffusion is more favorable for positive than for negative e. This is due to the 
fact that the used potentials are steeper in compression than in tension. 
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Figure 2.8: l c for (a) negative and (b) positive values of the misfit e. Only values of 
l c > 4 are taken into account. 



Finally we have a look at the mentioned critical island size l c and critical misfit e c 
from which on exchange is favored over hopping diffusion. Figure 2.8 shows the critical 
island size as a function of the misfit for the used pair-potentials. For both the positive 
and the negative branch of the misfit exchange diffusion becomes favorable at a critical 
misfit \e c \ which is the greater the steeper the potential is. At a given potential \e c \ 
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is greater for negative misfits than for positive misfits. This reflects the fact that the 
barriers for exchange are higher in this region of the misfit. 

For all considered potentials the critical island size starts to drop from \e c \ with in- 
creasing misfit In the high misfit region exchange becomes the determining diffusion 
mechanism already for small island sizes. 

In conclusion we find that at a given misfit a steeper potential leads to higher barriers 
for the exchange diffusion. This corresponds well to the fact that a steeper potentials 
delays the introduction of misfit dislocations [70] due to a higher activation barrier [19]. 

2.4 Calculations for the three-dimensional case 

In this section we present calculations for the 3d case. Due to the high computational 
demand and some conceptual problems regarding the application of the Molecular Static 
method to a 3d surface - which we will discuss later - the following results give a 
qualitative analysis of the different downward diffusion modes. We will also show that 
the 2d calculations indeed are of some relevance to the more realistic 3d case. 

Due to the isotropy of the used pair-potentials, particles interacting via equation 
(2.2) or equation (2.3) naturally arrange into a fee lattice structure. As interaction 
potential we choose the a = 5.0 Morse potential. Considering our 2d calculations this 
rather soft potential should display a strong dependency of the diffusion barriers on the 
misfit and the island size. 




Figure 2.9: Schematic representation of 
a hexagonal shaped island placed on a 
fcc(lll) surface (top view). Higher par- 
ticles appear in lighter grey. 



In our calculations we address the fcc(lll) surface. For symmetry reasons we place 
a hexagonally shaped adsorbate island on a substrate consisting of h = 11 layers (see 
fig. 2.9). One should notice that there are two different kinds of close-packed island 
edges on a fcc(lll) surface. Commonly the {100} microfacets of an island are labeled 
as A edges and the {111} microfacets are labeled as B edges. Due to the different local 
arrangement of the particles both kind of edges have to be treated separately for the 
calculation of diffusion barriers on descending steps. A particle is placed on top of the 
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island near the edge in question. In order to determine the minimum energy path for a 
descending move the test particle (for hopping: the particle on the top, for exchange: 
a particle from the island edge) is drawn in small steps parallel to the x-direction (see 
fig. 2.9). At each step the test particle is relaxed in the y- and z-direction. All other 
particles are relaxed without restrictions. A problematic point with 3d Molecular Static 
calculations is, that one can not be sure to really find the minimum energy path [71]. 
This is due to the huge numbers of possible particle positions compared to the 2d case. 
However, the method yields at least a qualitative analysis of the different diffusion 
mechanisms as a function of misfit and island size and allows a comparison of A or B 
steps. 

2.4.1 The homoepitaxial case 

We first analyze the barriers for descending diffusion moves from A and B steps in the 
homoepitaxial case (e = 0) for an island like shown in figure 2.9. Figure 2.10 shows the 
energy path for the four possible diffusion types: hopping and exchange diffusion on a 
A and B step, respectively. 



Hopping diffusion 

As figure 2.10 displays the energy paths for hopping diffusion down A and B steps look 
quite similar. Figure 2.11 (left panel) shows the way of a descending particle down an A 
step. The particle starts from a fcc-type binding site (minimum 1). Since the particle 
passes through an additional hep-type binding site on its way down the step a second 
minimum appears in the energy path (cf. fig. 2.10). At the end of the hopping move 
the particle is attached to the edge on a fcc-type binding site (minimum 3). 

The path of a particle downwards a B edge is similar (see fig. 2.11, right panel). The 
transition starts on a hep-type binding site, which is energetically slightly disadvanta- 
geous in comparison to a fcc-type site (minimum 1). A second minimum arises in the 
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Figure 2.10: Total energy of the 
system (substrate, island and test 
particle) for hopping and exchange 
diffusion from A and B steps in the 
homoepitaxial case (e — 0). 
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A-step 



B-step 



Figure 2.11: Left panel: 
Way of a particle down 
an A step. Right panel: 
Way of a particle down a 
B step. The numbers indi- 
cate positions of local en- 
ergy minima. Note that 
the particle starts from an 
fcc-type binding site on 
the A step and from a 
hep-type binding site on 
the B step. 



energy path, again due to a fcc-type binding site (minimum 2). As figure 2.10 shows 
particles attached to B steps are clearly weakly bound compared to particles attached 
to A steps. A similar result was found for different pair-potentials in the fcc(lll) ge- 
ometry before [62]. In conclusion - as a result from the similarity of the diffusion paths 
- the hopping barriers downward A and B steps are quite the same. 



Exchange diffusion 

The situation changes completely in the case of exchange diffusion. Figure 2.12 shows 
the paths for exchange diffusion on an A step (left panel) and a B step (right panel). 
For the exchange move on an A step the top of an underlying particle lies in front of 




A-step 




B-step 



Figure 2.12: Left panel: 
Exchange diffusion path 
on an A step. Right panel: 
Way of a particle down a 
B step. The numbers indi- 
cate positions of local en- 
ergy minima. 



the test particle (see fig. 2.12, left panel). This leads to the strongly enlarged barrier 
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in the energy path (fig. 2.10) for this diffusion move. 

On the B step the test particle faces a bridge site between two underlying substrate 
particles (see fig. 2.12, right panel). This leads to a rather small barrier for exchange 
diffusion on B steps, which is proved for various material types (see e.g. [62-68]). 



2.4.2 The heteroepitaxial case 

From the above made considerations it is clear that the arrangement of the particles in 
an A step is similar to the 2d situation: since the test particle has to overcome the top 
site of an underlying particle exchange diffusion is disadvantageous in the homoepitaxial 
case. Analogous to our findings for the 2d case the influence of the misfit (e ^ 0) is 
expected to be especially pronounced on exchange diffusion moves on A steps. Indeed, 
the energy path for e = 7% (see fig. 2.13, left panel) shows a strongly decreased 
exchange diffusion barrier for A steps. As figure 2.13 (right panel) in comparison to the 





X 

Figure 2.13: Left panel: Total energy of the system (substrate, island and test particle) 
for hopping and exchange diffusion from A and B steps in the heteroepitaxial case 
(e = 7%). Right panel: Way of a particle down an A step for e = 7%. 

homoepitaxial case (see fig. 2.12, left panel) displays this is due to the arrangement 
of the particles: the misfit causes a shift of the step particles over the top site of the 
underlying substrate particles. The arrangement of the A step particles becomes similar 
to that of B step particles, resulting in a rather low barrier for exchange diffusion. 



2.5 Conclusions and outlook 

In this chapter we investigated diffusion barriers from island edges. To this end a 
Molecular Static method was applied in two and three dimensions. The calculations 
have been performed for various Morse and the 12, 6 Lennard- Jones potential. Our 2d 
calculation have yielded the following results: 
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• The barrier for exchange diffusion moves down an island edge strongly depends on 
island size and misfit. The barrier decreases with increasing misfit and increasing 
island size, generally. 

• In the large misfit regime the barrier for exchange diffusion becomes smaller than 
the barrier for hopping diffusion. 

• The steeper the used pair-potential is, the larger is the critical misfit at which 
exchange diffusion becomes favorable. 

For the 3d case we performed calculations for the fcc(lll) surface using a Morse po- 
tential. Here, two results are of particular interest: 

• We have shown that the behavior of the diffusion barriers on A steps is qualita- 
tively the same as in the 2d case: an increasing misfit favors the exchange over 
hopping diffusion on A step. 

• In the large misfit regime exchange barriers on A steps become of approximately 
the same size as on B steps. 

Despite the fact that the Molecular Static method is a frequently used tool even in 
material specific barrier calculations, it implies some problems, especially regarding 
calculations in three dimensions. As mentioned above - considering the large number of 
particles one has to take into account - it is difficult to find the minimum energy path 
with this method. Given the interatomic potential the activation-relaxation technique 
(ART) [72-74] would perhaps be a more appropriate method for the calculation of 
diffusion barriers in three dimensions. 

One has also to bear in mind that both barrier and attempt frequency determine 
the rate for a diffusion event. Molecular Dynamics (MD) studies show that the attempt 
frequency for exchange diffusion is usually smaller than for hopping diffusion [75-79]. 
This favors hopping diffusion mainly at low temperatures. 
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Chapter 3 



Off-lattice Kinetic Monte Carlo 
simulations 

As described in chapter 1 atomic size mismatched systems always try to relax the 
strain imposed by the misfit. This can either be done by a rearrangement of the 
surface particles into mound-like structures or the introduction of misfit dislocations. 
Both relaxation mechanisms have in common that the distance between neighboring 
particles can vary over a wide range. For that reason the representation of the system 
by a rigid lattice is not longer practical for the simulation of heteroepitaxial growth. 
In order to model these relaxations in simulations it is crucial to allow for continuous 
particle positions. 

Given at least an approximation for the interatomic potential Molecular Dynamic 
(MD) simulations would be the most realistic way to do so (cf. chapter 1). Indeed there 
are several studies on the formation and mobility of misfit dislocations (see e.g. [48,51]). 
However, this method suffers generally from the restrictions to short physical time scales 
(< 10~ 6 s). In particular MBE relevant time scales of seconds up to minutes are even 
with modern computers not yet feasible. In conclusion, the most promising way to 
simulate heteroepitaxial growth today is by means of off-lattice Kinetic Monte Carlo 
(KMC) simulations. 

Historically, the first off-lattice KMC methods were the so-called ball and spring 
simulations, which are used with some success up to now [52-57]. In this kind of 
simulations strain is applied to the system by assuming a harmonic interaction (spring) 
between the particles (balls). However those models fail for example in the simulation of 
misfit dislocations and need somewhat artificial bond-counting rules for the calculation 
of binding energies (see chapter 1 for details). 

Here we follow a different direction based on work by Schindler and Wolf [61]. We 
apply the method to three different problems of heteroepitaxial growth. We assume 
that the interaction of the particles is given by a pair-potential which is a function 
of the continuous distance between the particles. The main idea of the method is to 
compute the activation barriers from the potential for each hopping event and to use 
the obtained rates in a standard rejection-free KMC simulation. In the following we 
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will discuss the method in detail along with its advantages and restrictions. 

3.1 Calculation of the activation energy 

We consider a system consisting of n particles interacting through a pair-potential 
Uij. The continuous distance between the two particles % and j is given by r^, where 
fij = (xij,Zij) T , fij = (xij,yij, Zij) T for the 2d, 3d case, respectively According to the 
symmetry implied by the potential the particles arrange into a crystal. Since all poten- 
tials used in this work (cf. appendix A) decrease rapidly towards zero with increasing 
particle distance a cut-off distance r cut with = for r\j > r cut is assumed. In most 
cases the cut-off is chosen to be r cut = 3r , where r is the equilibrium distance of the 
particles. 

As discussed in chapter 2 a particle deposited on the surface of a crystal moves in 
an energy landscape consisting of local minima (binding states) separated from each 
other by saddles (transition states) and maxima. One should notice that according to 
the definition of the potential energy surface (PES) given in chapter 1 in the 2d case 
the transition sites are given by maxima. However, if one considers the test particle's 
z-coordinate (perpendicular to the surface) as a free parameter the coordinates x and z 
along with the potential energy of the system reassemble again an energy surface where 
the transition site is given by a first order saddle point (see also fig. 3.1). In order to 
keep a clear notation in the following we will always refer to saddles when transition 
sites are considered. 

Since a particle hopping from one local minimum with energy E b to another has to 
overcome such saddles with energy E t , the goal is to calculate E b and E t and, hence, to 
obtain the activation energy E a = E t — E b for each diffusion event in a KMC simulation. 

3.1 .1 Searching for the saddle point 

Due to the influence of the adatom on the crystal one would actually need to relax the 
whole system in every step during the search for the transition site (like it was discussed 
in the last chapter). However, it was shown in [11,61] that a frozen crystal during the 
saddle search only shifts the activation barriers uniformly to slightly higher values 
(about 10% higher). Since calculations for a frozen crystal save a lot of computer time 
we restrict our method to this simplification during the calculation of E t . A further 
advantage is that during the barrier calculations one only has to keep track of the 
adatom energies, instead of the energy of the entire system. 

We now introduce the method used in our simulations for the calculation of the tran- 
sition energy. To this end an iterative algorithm for the saddle point search introduced 
in [72-74] is applied to our problem. This so-called activation-relaxation technique 
(ART) was originally developed to obtain energy minimized structures in glassy mate- 
rials: a system is allowed to evolve by following well-defined paths over saddles between 
local energy minima. From this method we implement the saddle search algorithm: 



32 



First the adatom is slightly displaced from its binding position in the wanted 
direction (i.e. the direction of diffusion). This causes the force F acting on the 
particle to become nonzero. 



By iterative application of the redefined force G given by 

G = F - (1 + a) (p£\ e 



(3.1) 



with a > the adatom is moved in small steps in direction of G toward the 
nearby saddle. Here e is the unit vector pointing from the last local minium to 
the current position of the particle. 

• The iteration ends when the saddle point is reached and G = F = 0. 

Since the redefined force G is opposite in sign to F in the direction parallel to e and 
equal to F in any direction perpendicular to e, the particle is forced in small steps a 
valley up-hill in the PES (see fig. 3.1). The positive number a controls the increment 
in each iteration step. As figure 3.1 shows, a too large value of a results in missing the 
relevant saddle. On the contrary, for small values of a, more iteration steps are needed. 
In our simulations a good value of a has to be found to reach the right saddle with as 
little iteration steps as possible. Values between a = 0.5 and a = 2.0 have yielded the 
best performance. 





- a - 


a=0.5 




a=4.0 


- B- - 


a=6.0 





Figure 3.1: Energy surface for the 2d case given as contour plot (left panel) and surface 
plot (right panel). The trace of the iteration (3.1) is shown for different values of a. 
Note that for a too high value of a (here for example a = 6.0) the method fails to end 
up in the relevant saddle. 

Formally, local maxima also fulfill the stopping criteria G = F = of the iteration 
(3.1). But the iteration can only end up in a local maximum if, by accident, a maximum 
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is chosen as starting point. Otherwise, the components of G in direction of F drive the 
adatom downhill away from the local maximum. 

In earlier simulations for the 2d different method has been used for the 

barrier calculation. As mentioned above for the 2d case the transition site is given 
- per definition of the PES - by a local maximum, with respect to the ^-coordinate 
only. Therefore E t can be calculated by minimizing the adatom's potential energy as 
a function of Zi at fixed horizontal coordinates Xi and maximizing this energy with 
variation of Xi. This is done by using a one-dimensional minimization method [69]. 

Both methods yield the transition site in comparable computing time. The iterative 
method has the advantage to work also for the 3d case. Moreover, the ART-like method 
proved itself more reliable in detecting the relevant saddle point. 

3.1.2 Calculation of the binding energy 

The calculation of the binding energy is computationally less demanding than the sad- 
dle search. In most cases it is sufficient to guess the position of the binding site by 
geometrical considerations. For example in the 2d case - when the particles arrange in 
a triangular lattice - at moderate misfit the atom in question will occupy the bridge site 
between two underlying particles. After the adatom is placed to the guessed position 
in a forthcoming relaxation of the system (see next section) the particle ends up at 
its binding site with E^. Only in the high misfit regime, when dislocations have to be 
considered this can fail. In this case, the particle is moved in direction of the binding 
site and a additional minimum search is performed. This can be done with the same 
methods as described for the saddle search. 

With the obtained binding and transition energies we are able to calculate the 
activation barrier for hopping diffusion rather fast. However, the described methods do 
not allow for the calculation of barriers of concerted moves, like the exchange diffusion 
examined in chapter 2. For this reason we restrict ourselves in the following to situations 
were these moves can be neglected. However, for a more comprehensive description of 
heteroepitaxial growth concerted moves will have to be incorporated in the simulation 
method in some way. 



3.2 Deformations of the crystal 

An adatom hopping on the crystal surface from binding site to binding site surely influ- 
ences the surface around its current position. Due to the additional binding partner the 
nearby particles will slightly change their positions. This effect is the more pronounced, 
the greater the strain in the system is. 

In order to account for these local deformations of the crystal, after each microscopic 
event - i.e. deposition and diffusion - the system is relaxed. This is done by minimizing 
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the total potential energy of the system, given by 



E to t = 



(3.2) 



, rij < r cut with respect to the particles' coordinates using a standard conjugate gra- 
dient method, taken from [69]. As a simplification, most of the time this relaxation 
is performed in a local way: only particles within a sphere of radius r cut around the 
locus of the last event are allowed to change their position during the procedure (see 
also fig. 3.2). This restriction saves a lot of computer time and is a valid simplification 
since the influence of a single microscopic event should be locally limited, due to the 
fast decreasing pair-potentials used in our simulations. 



However one may argue that especially at the edge of the considered sphere the 
local relaxation could add artificial strain to the system. For that reason a global 
relaxation with respect to all particle coordinates is performed after a certain number 
of microscopic events. This number depends on the specific problem (e.g. value of 
the misfit, temperature, deposition rates) and is guessed from preliminary simulation 
runs. As a rule of thumb the maximal change of the diffusion barriers due to the global 
relaxation should not exceed more than about 1%. 

One should notice that both mentioned types of relaxation do not lead to a sub- 
stantial rearrangement of the crystal but the coordinates and activation energies of the 
affected particles are changed slightly. The relative positions of the particles remain 
unchanged. 



3.3 Rejection-free Kinetic Monte Carlo simulations 



The binding energy and transition energy E tji have now to be calculated for each 
possible diffusion event i. The activation energy amounts to E a ^ = E t ^ — E^ (eq. 




Figure 3.2: Local relaxation in the 
2d circle of radius r cut is 

drawn around the active (here dark- 
est) particle. The (light gray) parti- 
cles within this circle are allowed to 
change their position during the pro- 
cedure. 
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(1.2)). Since we assume that the diffusion process follows an Arrhenius dynamics, see 

E a ^ 

chapter 1, the proper rates are given by Ri = u e~~i^. 

The obtained diffusion rates along with the rate for deposition of a new particle Rd 
are now used in a rejection-free KMC algorithm [4,80]. Because of the fact that espe- 
cially in the low temperature / high barriers regime the probability of non-modifying 
steps becomes significant non rejection-free techniques like the Metropolis algorithm 
result in a considerable slowing down of the simulation. Since in our simulations a 
diffusion barrier increases fast with increasing number of binding partners it pays to 
take the additional administration effort of a rejection-free continuous time algorithm: 
every iteration step of the algorithm a modifying event k is drawn and performed with 
the correct probability 

Then, the time interval r between two modifying steps is drawn randomly according to 
the probability of these steps. The time interval is given according to an exponential 
distribution by 

T = -\nt/(R d + J2 R *)> ( 3 - 4 ) 

i 

where £ is a uniformly distributed random number < £ < 1. Following the work of 
Ahr [4] the random selection of events is performed using a complete binary tree. 
Applied to our off-lattice method the algorithm works like this: 

1. An event k is drawn according to its probability pt (3.3) and performed. 

2. The crystal is locally relaxed around the location of the event. 

3. The rates for all diffusion events affected by the local relaxation (that is at least 
the particles within the relaxation sphere) are newly calculated. The search tree 
is updated accordingly. 

4. The system time is increased by the time interval r, equation (3.4). 

5. If the condition for a global relaxation is met the whole system is relaxed and all 
diffusion barriers are newly calculated. Accordingly the search tree is updated. 

6. The iteration starts from step 1 unless a stopping condition - for example when 
a certain system time is reached or a maximum number of particles has been 
deposited - is fulfilled. 

We should stress here that considering the used computer time the administration of 
events plays a minor role in our simulations. The time consuming part of our simulations 
is the calculation of the interaction potential between the particles. 
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3.4 Lattice based method 



Since calculations of the potential energy always involve a cut-off distance it is possible 
to predict at the beginning of each iteration run which particles will participate in 
the following energy calculations. So one could - when a event is drawn - search the 
whole crystal once for the relevant particles. In fact this saves a lot of computer time 
compared to a new online search through the whole crystal for every calculation of the 
interaction potential. But still all particles of the system have to be considered once 
per iteration which involves many expansive calculations of distances. 

However, in the case of moderate misfit, if the lattice structure is only bent but 
not broken, one can exploit the advantages of a lattice based method. This makes 
a double entry book-keeping necessary: on one hand each particle of the crystal is 
assigned to a lattice position in a perfect crystal - in the 3d case indexed by a triple 
of integer numbers (i,j,k). On the other hand each particle has still its coordinates 
given by real numbers (x, y, z). The condition for the adaptability of the method is that 
a non-ambiguous mapping between (i,j,k) and (x,y,z) remains possible throughout 
the whole simulation. That means in particular that the formation of dislocations is 
prohibited. 

The big advantage of the lattice based method is that the positions of the relevant 
particles are given by basic arithmetical operations of integer numbers. Calculations of 
the particles' real distances are only done when necessary which saves a lot of computer 
time during the simulations. A nice side effect is that one knows a lot about the 
structure of the growing crystal. This can, for example, be used to detect island edges 
during the simulations or allows to use pre-calculated diffusion barriers for recurrent 
situations. 

However, if one wants to simulate heteroepitaxial growth in the large misfit regime - 
depending on the choice of the potential - the lattice based method is inoperative. One 
way here to speed the simulations up is to divide the system into boxes of side length 
r cut . During the energy calculations for a certain particle only the box which contains 
this particle together with the neighboring boxes has to be considered. Thus only a 
small fraction of the system has to be searched for relevant particles. It is understood 
that this requires again a book-keeping over all boxes and their particles. This box 
method is applied in [70] to the off-lattice simulations and accelerates the simulations 
of dislocation formation significantly. 

3.5 Conclusions 

We have presented an off-lattice KMC algorithm which allows for the simulation of 
heteroepitaxial growth over a wide range of misfits. The following chapters will show 
that this simulation technique can be applied to various problems of heteroepitaxial 
growth and describes them quite successfully. 

We have to stress that within this method we are not yet aiming at the simulation 
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of distinct material systems but heteroepitaxial growth is examined in a qualitative 
way. To get a more realistic description of certain materials one would need to use 
material specific potentials. On the one hand such potentials for semiconductor (see 
e.g. [46,47]) or metallic (see e.g. [46]) systems are generally numerically much more 
demanding. Since this potentials are many body potentials and include long-range 
interactions (leading to a greater cut-off distance) they are currently out of scope. On 
the other hand, most of the material specific potentials are fitted to properties of the 
bulk like particle distances and elastic constants and most of them fail to describe 
properties of the crystal surface properly. 

In the following the aim is to gain general insights in strain-related phenomena 
rather than to obtain material specific results. We focus on phenomena observed in 
a various number of heteroepitaxial systems and therefore should not depend on a 
particular choice of potential. We restrict ourselves to fundamental questions like the 
influence of the misfit or the steepness of the used pair-potential on the described 
phenomena. 
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Chapter 4 

Simulation of misfit dislocations 



In this chapter we focus on the simulation of incoherent heteroepitaxial growth where 
misfit dislocations appear in the growing film. At moderate misfits <C 1) the 
adsorbate first grows coherently with the substrate i.e. the topology is that of a perfect 
crystal. However, the thicker the adsorbate film becomes, the higher becomes the elastic 
energy stored in the film. At the so-called critical thickness h d c the strain is relieved 
by the introduction of misfit dislocations (see also chapter 1). In this new incoherent 
state the crystal topology is perturbed near the substrate/adsorbate interface. With 
further deposition of adsorbate material the dislocations are buried and the lattice can 
grow with the adsorbate's lattice constant. In technical applications the formation of 
dislocations have desired as well as undesired effects. 

One example is the fabrication of II- VI semiconductor lasers which emit light in 
the blue and green region [81]. Here, for example ZnSe has to be deposited on a GaAs 
substrate implying a significant lattice mismatch of about 0.27% at room temperature. 
But the introduction of misfit dislocations can lead to unstable devices of poor reliability. 
In order to avoid the dislocation formation one needs to know the critical thickness under 
the given growth conditions and for the used alloy layers. 

On the other hand the impact of buried dislocation structures for example on the 
fabrication of magnetic nano-particles has been discussed, recently. One goal here is 
to grow patterned arrays of Co pillars on a metal substrate. To this end a thin Pt film 
is deposited on a sapphire substrate. The lattice mismatch between sapphire and Pt 
induces the introduction of dislocations in the Pt film. Since these dislocations repel 
each other the equilibrium configuration is given by a network of rather equally spaced 
dislocation arrays. If the Pt film is thin enough the variation of the lattice constant due 
to the dislocations causes modulations of the activation barriers for Co atoms diffusing 
on the film's surface. Since - like mentioned in chapter 1 - in metallic systems the 
diffusion is decelerated in regions of a tensile strained surface, the regular network of 
dislocations yields patterned nucleation of Co islands [82]. In this way the formation 
of dislocations opens a pathway for the self-assembly of novel magnetic data storage 
devices. 

In conclusion the knowledge about dislocation formation mechanisms, the depen- 
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dence on the properties of the used materials (especially the misfit) and the influence of 
formed dislocations on further growth of the crystal play an important role for technical 
applications. 

In this chapter we show that the introduced method (cf. chapter 3) is capable of 
simulating the formation of dislocations over a wide range of misfits. In the first part of 
the chapter we focus on the simulation of dislocation formation in the high misfit region 
|e| > 3%. We examine under which conditions the two mentioned types of dislocations 
- climb and glide (see chapter 1) - are formed. We determine the critical thickness as a 
function of the misfit and show that our results may be fitted well by a simple power law. 
In the second part of the chapter we address phenomena concerning buried dislocations 
for a region of relatively small misfit (1.4% < e < 2.2%). The simulation results on 
the evolution of the lattice constant as a function of the adsorbate film thickness are 
compared to experimental results, showing a good qualitative agreement. 



As described in chapter 3, the simulation of dislocations does not allow for the use of 
a lattice based method, which would reduce the computational effort a lot. For that 
reason we have to restrict our simulations in the following to rather small system sizes 
in 1 + 1 dimensions. Each simulation run starts with six atomic layers of substrate with 
a fixed bottom layer. Due to the fixed bottom layer no dislocations are introduced in 
the substrate during growth. The system size L (number of particles in the substrate's 
upper layer) is between L = 100 and L = 200. Within this range we found no significant 
dependence of the results on L. Adsorbate particles are deposited on the crystal's 
surface with a deposition rate Rj = lML/s. 

Because of its numerical feasibility all simulations are carried out for the Lennard- 
Jones 12, 6 potential 



(also see appendix A). We choose the same potential depth for all types of particle- 
interactions and set — U — 1.3125ey. In the homoepitaxial case (e = 0) this choice 
results in a rather high activation barrier for surface diffusion of about 0.90eV. In order 
to save computer time the interaction potential is cut off for distances > r cut 
with r cut = 3ro, where tq is the equilibrium distance between two nearest neighbors in 
the lattice. The interaction strength at r cut is less than 1% of the value at tq and can 
therefore be neglected. The interaction of two substrate particles is given by Uij (cr s ). 
Two adsorbate particles interact via Uy (a a ) whereas we assume that a substrate and 
an adsorbate particle interact via \ {U^ (<j s ) + (cr a ))- Measuring lengths in units of 
<r s , a a is chosen between 0.85 and 1.11, so we can simulate heteroepitaxial growth for 
misfits 



4.1 Growth in the high misfit region 




(4.1) 
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between —15% and +11%. For each value of e between 5 and 10 independent simulation 
runs are carried out. 

Our calculations in chapter 2 showed that in the case of the Lennard- Jones 12, 6 
potential one should take exchange diffusion for downward movement at step edges for 
values of the misfit e < —13% and e > 7% into account. However, it is conceptually 
rather complicated to implement concerted moves like exchange diffusion events in 
our off-lattice simulations. We therefore neglect exchange diffusion here. This is partly 
justified by the rather low temperature and high potential depths used in our simulations 
(T = 0.03^?), resulting in comparatively small rates for all interlayer diffusion processes. 

4.1.1 Formation of dislocations 

We now investigate the formation mechanism of dislocations. In our simulations we 
observe two different mechanisms of dislocation formation, leading to two different 
types of dislocations: climb and glide dislocations (see also chapter 1, fig. 1.4). 

Climb dislocations 

First we discuss the introduction of climb dislocations (see fig. 1.4(a)), characterized 
by a Burgers vector parallel to the substrate/adsorbate interface. To this end, figure 
4.1 displays the evolution of a crystal section for e = 10% - a rather high positive 
misfit, where the dislocations are found at the interface of substrate and adsorbate. 
The crystal section is shown for different coverage with adsorbate particles from XML 
(fig. 4.1 (T)) to 5 ML (fig. 4.1 (9)) in steps of half a monolayer. The vertical lines (A 
to D) in the samples indicate the final positions of the climb dislocations introduced 
during the growth. 

In the early stages of growth (fig. 4.1 0, (2)) the given section of the substrate 
is covered by three well separated adsorbate islands. The formation of islands is here 
due to the rather high Schwoebel barrier (cf. chapter 2) in combination with the low 
growth temperature. As one can see the islands fit the lattice spacing of the substrate 
rather well in the island's centers. Towards the rims of islands the adsorbate relaxes to 
its preferred lattice constant. This results in incoherent states at the rims indicating 
already the formation of dislocations. 

At a coverage of 2ML (fig. 4.1 (3), @) at point C two islands start to merge resulting 
in a climb dislocation at their contact point. In figure 4.1 (5) the same mechanism is 
observed at point B. Since the given samples are typical for the growth at high misfits we 
conclude that climb dislocations arise preferentially from the merging points of islands. 
Samples (B) to (9) show that in our model once formed climb dislocation do neither 
vanish nor move with further deposition of adsorbate particles. 

Figure 4.1 (9) shows the final state of the crystal at 5ML coverage. Note that 
the initial three islands can still be identified as the tops of mounds on the crystal's 
surface. Here the grey level for a particle indicates the particle's average distance 
to its nearest neighbors of the same kind: the lighter its grey level the more is this 
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Figure 4.1: Formation of climb dislocations: section of a crystal at a coverage with 
adsorbate particles between 1ML (T) and 5ML (9) in steps of half a monolayer. The 
adsorbate particles (e = 10%) appear in dark grey Final locations of dislocations are 
indicated by the vertical lines (A to D). In figure (9) the grey level of the particles 
indicates the particle's average distance to its nearest neighbors of the same kind. The 
lighter its grey level the more is this particle under compression. 
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particle under compression. In the case of e > adsorbate particles therefore appear 
darker near dislocations, since they are able to approach their preferred (greater) lattice 
constant there. In regions where the adsorbate grows coherently with the substrate the 
compressed adsorbate particles are displayed brighter. Note that the substrate is also 
influenced by the situation of the adsorbate film above: in case of positive misfit the 
substrate is under tension in regions of coherent growth, whereas near dislocations the 
substrate is compressed. 

Glide dislocations 




Figure 4.2: Formation of glide dislocations for a crystal sections, misfit is e = —5%. 
The coverage with adsorbate particles for the samples (T) to (8) is 12ML, 13ML, 15ML, 
16.75ML, YJML and 1BML, respectively. 



At lower misfits the resulting dislocations can be of a different kind. Here we often 
find the formation of glide dislocations - characterized by a Burgers vector with a 
component vertical to the substrate/adsorbate interface which does not contribute to 
the relaxation (also see fig. 1.4(b)). Figure 4.2 shows the evolution of a crystal section 
for a misfit e = —5% with a coverage of the substrate from Y2ML to 18ML. For the 
given value of the misfit dislocation formation does not start at the substrate/adsorbate 
interface but a few monolayers above. 
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Once again we observe that the formation of the dislocation starts in the valley 
between two mounds on the crystal's surface. Here the core of the dislocation is located 
- indicated by an arrow in sample (2). From the core the dislocation first grows in 
two branches shown in 4.2 (3), (4). With further growth (5) the left branch of the 
dislocation vanishes with only the right one remaining (5). Note that, although the 
dislocation is introduced a few monolayers above the substrate/adsorbate interface, 
substrate particles are influenced by the morphology of the adsorbate film: the substrate 
is slightly under compression in the region where the dislocations arise in the adsorbate 
film above. 

The reason for the misfit dependency of the formation mechanisms is due to the 
activation barrier for the introduction of a dislocation AEd- A recent study by Trushin 
et. al. [19] showes that in case of the Lennard- Jones system for misfits e > 8% the 
climb mechanism is the preferred dislocation formation mechanism. For smaller misfits 
both mechanisms - climb and glide - are competitive in energy costs. This findings are 
in good agreement with our simulation results, where in the large misfit regime climb 
dislocations are found to be the preferential dislocation type. 

As our simulations indicate dislocations of both types tend to form in valleys between 
two mounds respectively islands on the crystal surface. A rough surface therefore seems 
to facilitate the dislocation formation. Since particles in mounds can relax their strain 
up to a certain amount (see also chapter 5) the binding positions in the valleys between 
them are energetically especially unfavorable. For example in case of positive misfit 
particles in a mound's surface can move a little outwards resulting in a higher, more 
favorable distance to their binding partners. But binding sites in the valley where 
two mounds meet are then under compression and serve as seeds for the dislocation 
nucleation. Indeed simulations for smoother surfaces show [70] that the appearance of 
dislocations is shifted to higher thicknesses of the adsorbate film. 

4.1 .2 Number of dislocations 

In our simulations we observe that in each run several dislocations appear quite simulta- 
neously, within the range of a few monolayers. After the deposition of a few additional 
monolayers of adsorbate material the number of dislocations stays constant. It is clear 
that, the higher the value e of the misfit is, the more dislocations are needed in order 
to relax the strain stored in the crystal. In the following we discuss the functional 
dependence between the number of dislocations rid and e. 

To this end the number of dislocations is counted for each simulation run after 
the nucleation of dislocations has stopped. To prove the existence of a dislocation 
we determine the coordination number n c of each particle by calculating the Voronoy 
polyhedra [83,84]. Voronoy polyhedra are a generalization of the Wigner-Seitz cell to 
a system without a fixed lattice. The number of sides of a Voronoy polyhedron gives 
the coordination number n c (n c = 6 for a particle in a perfect triangular lattice). A 
Burgers circuit [16] is drawn around regions of the crystal with n c ^ 6. A non- vanishing 
Burgers vector then indicates the appearance of a dislocation. 
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Figure 4.3: Number of disloca- 
tions per system size rid/L as a 
function of the misfit e. The data 
points for the low misfit region 
(filled) are by courtesy of Ch. Vey 
The error bars represent as the 
standard error of the simulation re- 
sults. The dashed line gives the 
theoretical number of perfect dislo- 
cations in a system of size L. 



Figure 4.3 shows the number of dislocations per unit length no/L counted for each 
value of e. It is calculated after further deposition of six monolayers adsorbate material 
after the first appearance of a dislocation in the crystal. At this thickness of the ad- 
sorbate layer the maximum number of dislocations should be reached. The dashed line 
gives the theoretical number of dislocations in a system of size L under the assumption 
that no = L\e\ perfect climb dislocations have to appear on a rigid substrate in order 
to fit the adsorbate film onto the substrate. Perfect dislocations are those for which 
the crystal topology far from the substrate/adsorbate interface is the same as in the 
coherent state and the Burgers vector is therefore an integer multiple of the lattice 
vector. The formation of glide dislocations causes the deviations from the theoretical 
results for —7% < e < —3% and 4% < e < 8%. This is due to the fact that glide 
dislocations are spatially more extended than climb dislocations. For this reason in the 
case e > more and for e < less dislocations than n D = L\e\ have to be built when 
partial dislocations appear. The deviations for small values of the misfit e < 3% are 
mainly due to finite size effects and are discussed in detail in [70]. 

4.1.3 Critical thickness 

We now discuss the dependence of the adsorbate thickness h d c where dislocations first 
appear on the misfit as observed in our simulations. Since the high impact of disloca- 
tions on the quality of technical applications during the last 50 years several theoretical 
models have been proposed in order to calculate the critical layer thickness h d c for the 
formation of dislocations in heteroepitaxial growth. 

Energy balance model 

The first theoretical treatment was done by Frank and van der Merwe [85,86] in 1949. 
They compared the additional interfacial energy due to the formation of dislocations at 
the substrate/adsorbate interface with the elastic energy relieved by this dislocations. 
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The film thickness where both energies balance each other is then identified with h d . 
This calculation implies that only interfacial dislocations may appear. Since this ap- 
proach does only consider the energies of the relaxed and the strained state it neglects 
actual the mechanisms of dislocation formation. 

Force balance model 

In order to account for these mechanisms in 1974 Matthews and Blakeslee [17] proposed 
a force balance model which takes a dislocation nucleation mechanism of special impor- 
tance for technical applications into account: they assumed that preexisting dislocations 
in the substrate propagate into the adsorbate film (so-called threading dislocations). 
Once the force F £ , acting to elongate the threading dislocation in the interface due 
to the misfits, balances Fi, the dislocation line tension resisting the elongation of the 
dislocation, h d is reached [87]. This yields an implicit equation for the critical thickness 



where b is the Burger's vector and C is a constant depending on the orientation of the 
dislocation line and the elastic properties of the adsorbate material. 

Both the thermodynamic and the force balance approach yield identical results for 
the calculation of h d [18,48] as a function of material parameters like the shear modulus 
or the Poisson ratio. However, since for a misfit greater than 4% the critical thickness 
becomes of the order of the lattice constant these continuous approaches are no longer 
valid and the discrete nature of the atomic layers should be taken into account [15]. 
Furthermore even in the low misfit regime experimental values for h d c are usually larger 
then the calculated ones [15,48,81,88]. Possible reasons for this poor agreement (with 
discrepancies as high as an order of magnitude) are that the continuous approaches take 
only dislocations near the substrate/adsorbate interface into account, different growth 
conditions (e.g. varying flux or temperature) are disregarded and - most importantly 
- kinetic barriers for the introduction of dislocations are neglected. As was shown 
in [19] the incoherent, relaxed state is separated from the coherent, strained one by 
a activation barrier AEj for the introduction of dislocations. This barrier stays finite 
even if the relaxed state is energetically favorable. The experimentally observed h d 
may therefore exceed the calculated values. Experimental results on the temperature 
dependence of the critical thickness also support the idea of strain relaxation as an 
activated process [89,90]. 

A power law for the critical thickness 

In order to compare our results on the critical thickness with experimental works we 
follow here a different approach proposed by Cohen-Solal and coworkers [88,91]. There, 
an energy balance model is proposed for calculating the critical layer thickness in het- 
eroepitaxial growth of semiconductor compounds. To this end the classical strain en- 
ergy, without any change of the substrate or dislocation formation, and the deformation 
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energy due to a full system of interfacial misfit dislocations were compared. The ener- 
gies were calculated using Keating's valence force field approximation [88]. The method 
yields a power law dependency of h d on e: 

h d = a*e~ 3/2 , (4.4) 

where a* is a material specific fit parameter. 

Despite the fact that the model neglects both the mechanism of dislocation nucle- 
ation and the kinetics of strain relief it shows excellent agreement with experimental 
data. The authors themselves fitted the power law to measured critical thicknesses for 
IV-IV, III-V and II- VI semiconductor compounds revealing an excellent agreement 
with values of a* between a* = 0.12 and a* = 0.50. Pinardi et al. [81] demonstrated 
in an independent study that equation (4.4) shows much better agreement with ex- 
perimental data for II- VI semiconductor compounds (with a* = 0.45) then the above 
described force balance method. 

Simulation results and discussion 

Figure 4.4 shows the critical layer thickness h d plotted versus the absolute value of the 
misfit e. For —0.03 < e < 0.02 the critical thickness is too large to be observed in our 



Figure 4.4: Critical thickness h d c 
versus misfit \e\ for e < (upper 
curve) and e > (lower curve). 
The error bars are obtained as the 
standard error of the simulation re- 
sults. The solid lines are calculated 
using equation (4.4) with a* = 0.15 
for e < and a* = 0.05 for e > 0. 
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simulations. The simulation results show a strong dependence of h d on the sign of the 
misfit. This was found before by L. Dong et al. [48]. We believe this dependence is due 
to the fact that, e.g., the Lennard- Jones potential is not harmonic. The potential is 
steeper in compression (e > 0) than in tension (e < 0), so that for e > it becomes 
favorable to form a dislocation for smaller values of h d . This corresponds well to the 
fact that in case of the 12, 6 Lennard- Jones potential the activation barrier for the 
introduction of dislocations AE d is found to be generally higher for adsorbate films 
under tension then for those under compression [19]. 
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Our simulation results agree well with the power law equation (4.4) (see solid lines 
in fig. 4.4). A nonlinear fit of our results yields a* = 0.15 for e < and a* = 0.05 for 
e > 0. 

Thus, our quite simple 1 + 1 dimensional model shows qualitatively the same de- 
pendence of the critical layer thickness on the misfit as semiconductor samples grown 
in molecular beam epitaxy. The values of a* obtained from our simulations show the 
same order of magnitude as the experimental ones but are expected to depend on the 
dimension and the applied interacting potential [70,88,91]. 

4.2 Simulations for the low misfit region 

In this section we examine the development of the lattice constant during heteroepi- 
taxial growth for the region of low positive misfits (e < 3%). Of special interest is the 
dependency of the lattice constant on the film thickness and the misfit. The examina- 
tions are mainly motivated by experimental work for the strained layer growth of ZnSe 
on GaAs. 

4.2.1 Development of the lattice constant during ZnSe/GaAs het- 
eroepitaxy 

Recently it has become possible [92, 93] to monitor the vertical lattice constant a± 
(lattice constant in the direction of growth) averaged over the whole adsorbate film 
during MBE growth. To this end a novel X-ray diffraction method was used which 
allows - unlike conventional methods - for the real-time in-situ X-ray diffraction (RIX). 

Conventional methods for the X-ray diffraction - like the u—29 scan at a symmetrical 
Bragg reflection - have two important disadvantages which restrict their application 
during MBE growth. First an extremely exact sample adjustment is required which 
is barely achievable under growth conditions. More importantly, a time-consuming 
angular scan by rotating both sample and detector has to be performed which takes up 
to 5 minutes per shoot, which simply eliminates real-time monitoring. 

The new RIX method circumvents these problems by using an asymmetric Bragg 
reflection (here the 113 reflection, resulting in an especially small exit angle) which 
yields the same results for a 26 scan of the reflection. In order to avoid the movement 
of the sample/detector during the measurement a slightly divergent X-ray beam is 
used which is detected by a multichannel CCD camera (also see fig.4.5 left panel). This 
enables to cover the whole 28 range within a few seconds by taking a single image with 
the CCD camera. The method yields two peak positions for the adsorbate and the 
substrate, respectively. 

By analyzing the distances between both peaks the vertical adsorbate lattice con- 
stant, averaged over the whole adsorbate film can be determined. Together with the 
growth rate, evaluated by RHEED (Reflection High Energy Electron Diffraction) oscil- 
lations, this gives a± as a function of the film thickness. 
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Figure 4.5: Left panel: Experimental setup for RIX, consisting of a conventional 
MBE chamber with a standard X-ray source, providing slightly divergent X-rays and 
a CCD camera. Right panel: Development of the ZnSe film averaged vertical lattice 
constant with increasing layer thickness as obtained from RIX. The curve shows the 
initial coherent growth and the subsequent gradual progress of relaxation. Data is 
provided by courtesy of A. Bader [92]. 



Figure 4.5 (right panel) displays the development of dj_ as measured for the growth 
of ZnSe on GaAs for a growth temperature of 520K. It can be seen that in the early 
stages of growth a± first remains approximately constant. This is believed to be due to 
pseudomorphic (coherent) growth of substrate and adsorbate in this region. The lateral 
compression of the larger ZnSe adsorbate (the misfit is e = 0.31% at growth tempera- 
ture) causes a increased vertical lattice constant of 0.57nm which is in agreement with 
this misfit. At a layer thickness h of about 200ML (llOnm) a± begins to drop. This 
indicates the relaxation of the ZnSe film due to the introduction of misfit dislocations: 
due to the reduced lateral compression of the film the vertical lattice constant moves 
towards the smaller bulk value. 

Note that the introduction of dislocations is guessed here by the development of a± 
and can not be observed directly during growth. In order to gain a deeper insight in 
the relaxation mechanisms leading to the observed behavior of a± we now adapt our 
simulation method to meet the requirements of the described problem. 

4.2.2 Adaptation of the simulation method 

We have to realize an expanded region of coherent growth in our simulations. As 
described above this is in principle possible by choosing lower values of |e|, yielding a 
higher critical thickness. However, misfits of the order of magnitude like those found 
for ZnSe/GaAs growth would require very large system sizes in order to avoid finite size 
effects in our simulations. Since the computation time necessary for these simulations 
would be beyond our present capacity we choose a different route here. 
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On the one hand the misfit is chosen moderately low (1.4% < e < 2.2%) requiring 
a lateral system size of L = 400, which is well within our computational range. On the 
other hand the surface is smoothed by applying a downward funneling method [4, 94, 
95] to the algorithm. Since, like mentioned, a rough surface promotes the formation 
of dislocation this helps to produce thicker coherent films without causing a severe 
deceleration of the calculations. 




Figure 4.6: Schematic representa- 
tion of downhill funneling: a newly 
deposited particle (dark) jumps suc- 
cessively to lower neighboring binding 
sites. In case of two equal low sites at 
the beginning of the process the direc- 
tion of the first jump is randomly cho- 
sen. 



Figure 4.6 shows a newly arrived particle jumping successively to lower neighboring 
binding sites until the lowest one is reached. The downhill motion is justified by the 
assumption that a newly deposited particle has a considerably higher kinetic energy. 
This leads to an increased mobility until the excess energy is dissipated. At the end 
of this process the particle gets bound at an energetically favorable site with a high 
coordination number. Results of Molecular Dynamics simulations (cf. chapter 1) con- 
firm this considerations [96,97]. Various growth models showed that an incorporation 
of the funneling mechanism yields correct results for the topology of growing surfaces 
(see e.g. [4,5]). 

With respect to the higher film thicknesses the substrate thickness is increased to 
11 layers in the following. All other settings of the algorithm remain unchanged. Ch. 
Vey [70] realized the thus modified algorithm and performed the simulations. 

In the following we present results for misfits 1.4% < e < 2.2%. For each value of 
e the results are averaged over 10 independent simulation runs. First we present the 
development of (2j_ SIS db function of the film thickness and discuss the different stages of 
growth like observed for one value of e. Subsequently we examine the influence of the 
misfit on the averaged vertical lattice constant. 

4.2.3 Development of the lattice constant with the film thickness 

Figure 4.7 shows the development of a± with increasing film thickness for a misfit 
e = 1.6%. The characteristics of the curve are in qualitative agreement with the 
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experimental data (see fig. 4.5, right panel): first the averaged vertical lattice constant 
remains constant followed by a subsequent decline of a±. 
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Figure 4.7: Development of the 
averaged vertical lattice constant 
with increasing layer thickness for 
e = 1.6%, showing the initial co- 
herent growth (a), the subsequent 
gradual progress of relaxation (b) 
and further growth of the relaxed 
adsorbate film (c). The dashed line 
gives the value of a c ° nst according 
to equation (4.5). The error bars 
represent the standard error of the 
simulation results. 



In case of e = 1.6% the constant region persists for about 80ML. Figure 4.9(a) 
shows that within this film thickness the adsorbate indeed grows pseudomorphic with 
the substrate without the introduction of dislocations. Due to the compression of the 
adsorbate film in lateral direction the vertical lattice constant is shifted to higher values 
as one would expect for the bulk crystal. The measured value of a± can be calculated 
approximately by a simple geometric approach. Consider two particles in a closed layer 



of the compressed film and a third particle located in the binding position between them 
(fig. 4.8). Due to the pseudomorphic growth the lateral distance of the two particles 
within the layer is given by r - the equilibrium distance of substrate particles. The 
particle on top is free to choose the preferred distance r (l + e) from its underlying 
nearest neighbors. The distance A between the uppermost particle and the underlying 
layer results then to 



Since all layers in the coherent region of growth should have approximately the distance 
A from each other this gives a estimated value for d± as function of the misfit in this 
region. Indeed the thus calculated a c f nst shows a good agreement with our simulation 
results (dashed line in fig. 4.7). 

As our simulation results verify (see e.g. fig. 4.9(b)) the region of decreasing a± is 
governed by the introduction of misfit dislocations. In the coherent film the formation of 




Figure 4.8: Calculation of the dis- 
tance A between a particle and the un- 
derlying compressed layer. 




(4.5) 
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(a) (b) (c) 

Figure 4.9: Section of a simulation run for e = 1.6% (the substrate is not shown). The 
three stages of growth: (a) coherent growth {h = 70ML), (b) formation of dislocations 
(h = 85ML) and further growth of the relaxed adsorbate film (h = 200ML). 

dislocations leads to a reduction of the compression and therefore to higher distances of 
the particles within the film. According to our simple model this decreases the distance 
between subsequent layers and yields thus the reduction of the averaged vertical lattice 
constant. 

In the third stage of growth after the strain in the adsorbate film is relaxed and 
the introduction of dislocations has finished (see fig. 4.9(c)) the distance between two 
subsequent layers A is that of bulk adsorbate. Since a± is averaged over the whole film 
(compressed and relaxed parts) a± slowly approaches the adsorbate bulk value. 

In conclusion we are able to reproduce the experimental results on d± as a function 
of the film thickness qualitatively. Our results confirm the explanation given in [92] for 
the development of a±_ with the thickness of the adsorbate film. 



4.2.4 Dependence of the lattice constant on the misfit 

We now focus on the dependence of a on the misfit between the substrate and the 
adsorbate layer. To this end we determine the value of a c l nst for 1.0 < e < 2.4. As 
figure 4.10 shows a c ° nst depends linearly on the misfit. Considering equation (4.5) this 
is easy to understand. Since we are dealing with rather small misfits here equation (4.5) 



52 



I CO 



0.996 ~ 
0.994 - 
0.992 - 
0.99 
0.988 - 
0.986 - 

I 

0.984 - 
0.982 - 
0.98 

0.978 - 

c; ' 
0.976 - 



J3 



0.01 0.012 0.014 



0.016 0.018 

e 



0.02 0.022 0.024 



Figure 4.10: Development of the 
averaged vertical lattice constant 
function of e. The dashed line 
gives the values of d± according to 
equation (4.6). 



can be linear approximated with 




As mentioned above in the coherent film A should be the same for all layers. Equation 
4.6 gives therefore an approximation for the functional dependence of a c f nst on e. The 
dashed line in fig. 4.10 shows the calculated vertical lattice constant which agrees well 
with the simulation data. Figure 4.11(a) shows the development of d± scaled with a c f nst 
for 1.4 < e < 2.2. As one would expect the region of decreasing a± is shifted to higher 
values of the film thickness the smaller the misfit becomes. Since this stage of growth 




Figure 4.11: Development of a±/d c f nst with increasing layer thickness for different 
values of the misfit. The right panel (b) shows the same curves with the film thickness 
scaled by s~ 3 / 2 . 

is governed by the critical thickness h d c for the formation of dislocations it is obvious to 
assume that the showed curves scale according to equation (4.4) with e~ 3 / 2 . As figure 
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4.11(b) shows this guess leads to a rather good collapse of the curves in this region of 
the film thickness. The deviation for e = 1.4% is believed to be due to finite size effects 
(see [70] for details). 

4.3 Conclusions and outlook 

In conclusion we have shown that within the framework proposed in chapter 3 it becomes 
possible to examine the appearance of misfit dislocations as one possible relaxation 
mechanism in heteroepitaxial growth by the means of KMC simulations. 

It was found that - like experimentally observed and theoretically predicted - also 
in our simulation dislocations appear at a certain thickness of the adsorbate film h d c . 
We detected that h d c is linked to the misfit by a simple power law, which also was 
measured for various semiconductor compounds. Despite the simplicity of the model 
(1 + 1 dimensions, simple pair-potential as particle interaction) it was possible to gain 
deeper insight into dislocation formation mechanisms. Furthermore we were able to give 
reasons for the experimentally observed development of the vertical lattice constant as 
a function of the film thickness. 

But we also encounter technical problems which limit the application of the algo- 
rithm on a realistic treatment of e.g. the dislocation evolution. Like mentioned above 
our model is not capable of treating concerted moves. Therefore it does not allow for 
the calculation of the activation barriers of dislocation movement, which is often ob- 
served in real materials. Whether it is conceptually possible to solve this problem is 
not yet clear. 

However, there are numerous further problems concerning dislocations in heteroepi- 
taxial growth which can be addressed within the algorithm. One important question 
is certainly the influence of buried dislocations on the topology of the crystal surface 
far from the disturbance. Preliminary results show that over-grown dislocations influ- 
ence the evolution of the adsorbate over a wide range of the film thickness. Another 
interesting problem is the influence of the steepness of the used potential on the dis- 
location formation. Especially the dependence of h d c on the used potential should be 
clarified. Since real MBE growth mostly takes place on stepped substrates this should 
be also examined in simulation which is quite easy to realize in our algorithm. One 
important mechanism for the the misfit dislocation formation in real materials is the 
propagation of preexisting substrate dislocations into the adsorbate film. The influence 
of those threading dislocations on the critical thickness is also a problem which could 
be addressed by the method. 

Last but not least the algorithm should be extended to 2 + 1 dimensions. Here 
one will encounter the problem of finding the relevant transition states for the diffusion 
jumps. It is also not clear if it is possible to reach the relevant system size and time 
scales for this computationally demanding problem. However, only simulations in 2 + 1 
dimensions will open the pathway to more realistic lattice structures, realistic potentials 
and thus more material specific predictions. 
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Chapter 5 



Simulation of 

Stranski-Krastanov-like growth 

As mentioned in chapter 1 dislocations (cf. chapter 4) clearly dominate the strain 
relaxation in sufficiently thick films and for large misfits. In material combinations 
with relatively small misfits an alternative effect can govern the initial film growth: if 
the system manages to overcome the barrier E 2 d-3d instead of growing layer by layer the 
adsorbate aggregates in 3d islands. The term island is commonly used to indicate that 
these structures are separated - unlike the emergence of mounds due to the Ehrlich- 
Schwoebel (ES) or the Asaro-Tiller-Grinfeld (see e.g. [15]) instability. 

In chapter 1 we already discussed the two possible growth modes, known for the 
separation of adsorbate material into islands: the Volmer-Weber growth mode with 
islands forming on uncovered substrate and the Stranski-Krastanov (SK) growth mode 
where 3d islands are found upon a pseudomorphic (dislocation free) wetting-layer of 
adsorbate material. This SK growth is observed for various strained heteroepitaxial 
systems, always for a quite large, positive misfit (2% < e < 7%). 

In order to avoid conflicts with other definitions and interpretations of the SK growth 
mode in the literature we will resort to the following specifications: 

• The adsorbate film growth in a layer-by-layer way pseudomorphic up to a kinet- 
ically controlled wetting-layer thickness h*. 

• 3d islands appear suddenly and quite simultaneously, marking the so-called 2d— 3d 
transition. 

• The further growth of the 3d islands is fed both by capture of newly deposited 
adsorbate particles and the decomposition of the supercritical wetting-layer. 

• Growth results in well separated 3d islands of similar shape and size, on top of 
the stable wetting-layer with reduced stationary thickness h c of a few monolayers 
height (h c < 7 ML). 
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Besides this basic processes a variety of phenomena play important roles for SK growth 
in real systems. For example the interdiffusion of adsorbate and substrate material or 
the segregation of compound adsorbate are believed to be of high relevance in many 
cases [98,99]. However, SK-like growth is observed in a variety of material systems 
which may or may not display this specific features. For instance, intermixing or seg- 
regation should be irrelevant in the somewhat exotic case of large organic molecules 
like PTCDA deposited on a metal substrate, e.g. Ag(lll). Nevertheless, this sys- 
tem displays SK growth in excellent accordance with the above specifications [100]. 
Most prominent examples for SK systems are IV-IV (see e.g. [101-107]), III-V (see 
e.g. [13,108-120]) and II- VI (see e.g. [121-123]) semiconductor compounds, see e.g. for 
InAs/GaAs self-assembled quantum dots figure 5.1. Despite the extensive investigation 



of SK growth, a complete detailed theoretical picture is still lacking, apparently. This 
concerns in particular the nature of the 2d — 3d transition. One problem clearly lies 
in the richness of the phenomenon. On the other hand, the mere diversity of SK-like 
systems gives rise to the hope that this growth scenario might be governed by a few 
basic universal mechanisms. Accordingly, it should be possible to capture and identify 
these essential features in relatively simple prototype systems. 

This hope motivates the investigation of simplifying models without aiming at the 
reproduction of material specific details. Some of the key questions in this context 
are: under which conditions does a wetting-layer emerge and persist? How does its 
thickness before and after the 2d — 3d transition depend on the growth conditions? 
Which microscopic processes trigger and control the sudden formation of 3d islands? 
How do the island size and their spatial arrangement depend on the parameters of the 
system? 




Figure 5.1: 200 x 200 nm scan- 
ning tunneling microscopy image of 
InAs/GaAs self-assembled quantum 
dots [124]. 
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5.1 Simulation model 



The method introduced in chapter 2 is now adapted for the the simulation of SK-like 
growth in 1 + 1 dimensions. The pairwise interactions between particles of the system 
are given according to a 12, 6 Lennard-Jones potential (see also appendix A for details) 



The choice of the parameters E^cr^- in equation (5.1) characterizes the different ma- 
terial properties in our model: interactions between substrate (adsorbate) particles are 
specified by the sets U s ,a s and U a ,a a , respectively. In principle, an independent pair 
of quantities U as ,a as could specify the substrate-adsorbate potential. For the sake of 
simplicity, aiming at a low number of free parameters, we use 



and a as = (a s + a a ) /2 for the interaction between the two species of particles. 

As the equilibrium distance r between two Lennard-Jones particles in the bulk is 
proportional to <Jij the lattice misfit e becomes e = (a a — a s ) ja s . Since SK growth 
in real systems is almost exclusively observed for positive misfits [27], we consider 
only cases with a a > a s here. If not otherwise specified, we have set the misfit to 
e = 4% - a typical value for SK systems like, e.g., in Ge/Si heteroepitaxial growth. As 
demonstrated in chapter 2 we have to take only hopping diffusion into account for such 
a small misfit. Exchange diffusion events can be neglected. In order to save computer 
time the potential Uij is cut off at a distance > 3r where the interaction strength 
is less than 1% of the value at the equilibrium distance. 

In our model growth takes place on six atomic layers of substrate with a fixed 
bottom layer and periodic boundary conditions in horizontal direction. The system size 
L (number of particles per substrate layer) is set to L — 800 in the following. Earlier 
simulations of smaller system sizes (L = 400, 600) revealed no significant L-dependence 
of the results presented here. 

As we have demonstrated in chapter 4, strain relaxation through dislocations is not 
expected for the chosen misfit within the first few monolayers. Since we deposit in the 
following only a few monolayers of adsorbate per simulation run, we are able to apply 
the latticed based method (chapter 3). 

An important modification concerns interlayer diffusion. Lennard-Jones systems in 
1 + 1 dimensions display a strong Ehrlich-Schwoebel barrier which hinders such moves 
at l<i-island edges (see chapter 2). This barrier is by far less pronounced in 2 + 1 
dimensional systems, because interlayer moves follow a path through an energy saddle 
point rather than the pronounced maximum at the island edge. In our investigations 
of the SK scenario we remove the Schwoebel barrier for all interlayer diffusion events 
by hand. One motivation is the above mentioned over-estimation in 1 + 1 dimensional 
systems. More importantly, we wish to investigate strain induced island formation 




(5.1) 



u as = ^ujr a 



(5.2) 
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without interference of the ES instability. Note that the latter causes the formation of 
mounds even in homoepitaxy [3, 15]. 

5.2 Stable wetting-layer thickness 

We will now first investigate the stable wetting-layer thickness h c , which is adjusted 
after the process of 2d — 3d transition is finished. For semiconductor systems wetting- 
layer thicknesses from h c ~ XML for the InAs/GaAs system [119] up to h c ~ AML 
(e.g. Ge/Si [103]) or even h c p=j QML for GalnAs/GaAs compounds [120] have been 
measured. Sometimes much thicker wetting-layers (h c > 50ML) are reported - however, 
with respect to the above given definition we do not refer to such cases as SK growth. 

A kinetic reason for the formation of a stable wetting-layer is the fact that the 
diffusion of adsorbate particles is found to be slower on the substrate then on subsequent 
adsorbate layers. For example in case of Ge growing on Si [103] the activation energy 
for the diffusion of Ge on the substrate is about O.QAeV, whereas the barrier for Ge 
diffusing on the first adsorbate layer is about OAOeV. A similar behavior is reported 
for the InAs/GaAs heteroepitaxy. Here the activation energy for the desorption of In 
located on the top of the wetting-layer is determined as 3.2eV and direct desorption 
from the substrate corresponds to a barrier of 3AeV [119]. Again the adsorbate particles 
are better bound to the substrate than to a pseudomorphic strained layer of the own 
species. In conclusion - at least in some semiconductor systems - the wetting-layer 
is stabilized by a slower diffusion on the substrate and a faster diffusion on top of the 
wetting-layer, as deposited particles will reach and fill gaps on the substrate more easily. 

These findings are also confirmed by our simulation results. For U as < U a the 
adsorbate layer thickness decreases from h* to h c = with adsorbate islands located 
directly on the substrate, eventually. For U as > U a a stable wetting-layer thickness 
h c > persists. If we set - at a simulation temperature T = 500i^ - U a — 0.74eV, 
for instance, the choice U as = 2.7eV (i.e. U s = lOeV, according to eq. (5.2)) leads 
to a layer thickness h c ps 2ML, whereas for U as = 0.8QeV (i.e. U s = l.OeV) a single 
monolayer of adsorbate forms the stationary wetting-layer (also see fig. 5.2). 

Figure 5.3 shows the corresponding activation energies for the diffusion of an adsor- 
bate particle on the substrate and the two subsequent adsorbate layers as a function 
of U as . As one can see in the case U as = 0.86eV the barrier for diffusion on the first 
wetting-layer becomes E a 0A7eV, whereas the barrier for diffusion on the substrate 
is given by E a ps 0.57eV. For U as = 2.7eV similar values are found for the diffusion 
on the second and first adsorbate layer, respectively. We find that an increase of the 
activation barrier for diffusion on two subsequent layers by roughly O.leV results in 
the formation of one additional stable adsorbate layer. In principle this behavior could 
extend to more then just two adsorbate layers. However, due to the short-range nature 
of the Lennard- Jones potential the influence of the substrate (given by eq. (5.2)) es- 
sentially vanishes on wetting-layers of three or more layers thickness. In the following 
we will use U as = 0.86eV (i.e. U s = l.OeV) and restrict ourselves to h c ps 1ML as a 
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Figure 5.2: Snapshots form 
simulation results for U as = 
0.86eV (upper panel) with 
h c « XML and U as = 2.72eV 
(lower panel) with h c ~ 2ML, 
both at T=500K. The six bot- 
tom layers correspond to the 
substrate. The darker the grey 
level of a particle, the big- 
ger the average distance to its 
nearest neighbors of the same 
particle type. 
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Figure 5.3: Diffusion barrier for ad- 
sorbate diffusion on subsequent lay- 
ers on the surface function 
of the adsorbate-substrate interaction 
strength U as . 



prototype SK scenario. However, the results presented in this work should not depend 
on the precise value of h c , qualitatively. 



5.3 Diffusion process 

Before analyzing the SK-like growth scenario, we compare the barriers for hopping 
diffusion in various settings on the surface. One should note here that on the one 
hand investigations of systems like Ge/Si(001) reveal a very complicated scenario due 
to anisotropies and the influence of surface reconstructions [13, 103, 125] which can 
not be modeled within our 1 + 1 dimensional system. On the other hand for most 
semiconductor systems (e.g. Ge/Ge(lll) [12] or In/GaAs(001) [13]) the barrier for 
hopping diffusion is higher on the surface of a compressed crystal, whereas diffusion is 
faster on the relaxed crystal. 
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5.3.1 Influence of the island height on diffusion barriers 



As mentioned in chapter 1 pair-potential systems dot not reproduce this feature in gen- 
eral, since mechanical compression of the crystal lowers the barrier for surface diffusion 
(see also fig. 5.5). However, in the mismatched two species system, it is more important 
to compare diffusion on the substrate, the wetting-layer and the surface of partially re- 
laxed islands. As shown above, the slow diffusion on the substrate in comparison to the 
fast diffusion on top of the wetting-layer stabilizes the wetting-layer. 
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Figure 5.4: Diffusion barriers as ob- 
tained in our model for a single adatom 
on a flat symmetric multilayer island 
with 24 base particles and one to six 
layers height. The symbols represent 
the activation energies for hops from 
the particle position x to the left neigh- 
bor site. The island is placed on top of a 
single wetting-layer. The leftmost bar- 
riers correspond to downward jumps at 
the island edge with suppression of the 
Schwoebel effect. Horizontal lines mark 
the barriers for adatom diffusion on the 
wetting-layer and on fully relaxed ad- 
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For the SK-like scenario the diffusion on islands of finite extension is particularly 
relevant. Figure 5.4 shows the barriers for diffusion hops on islands of various heights 
located upon a wetting monolayer. Note two different features: 

• Diffusion on top of islands is, in general, slower than on the wetting-layer and the 
difference increases with the island height. In our model, this is an effect of the 
partial relaxation or over-relaxation in the top layer of the island. 

• Depending on the lateral island size and its height, there is a more or less pro- 
nounced diffusion bias towards the island center, reflecting the spatially inhomoge- 
neous relaxation (cf. chapter 1). Note that this effect has to be distinguished from 
the diffusion bias imposed by the Schwoebel barrier, which would be present even 
in homoepitaxy and with particle positions restricted to a perfect undisturbed 
lattice. 

Clearly both features favor the formation of islands upon islands and hence play an 
important role in our simulation of SK-like growth. They concern adatoms which 
are deposited directly onto the islands as well as particles that hop upward at edges, 
potentially. As we will argue in the following section, upward diffusion moves play the 
more important role for the 2d — 3d transition. 
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5.3.2 Influence of the misfit on diffusion barriers 



0.51 1 1 1 1 1 1 1 

0.505 -X 

0.495- 

> 

g 0.485- 

■i- aJ v © 

LLI o.48 



0.475 
0.47 

0.465 
0.46- 



Figure 5.5: Diffusion barriers E]; 
for an adatom on the first wetting-layer 
decreases linearly with increasing misfit 
e. 



0.01 0.02 0.03 0.04 0.05 



Besides the island height also the misfit is believed to have a strong influence on the 
diffusion behavior of an adatom diffusing on an island. For that reason we measure the 
diffusion barrier for a single adatom on a flat symmetric island, located on one layer of 
adsorbate, for misfits 0% < e < 4.5%. Note, that with increasing misfit diffusion on the 
first wetting-layer becomes faster, since the adatom feels a more compressed adsorbate 
layer (cf. chapter 1). Figure 5.5 shows the diffusion barrier E\ WL for a particle on 
the first wetting-layer as a function of the misfit. As mentioned above the activation 
energy decreases linearly with increasing misfit (also see [11]). 



- ► 


4.5% 




4% 




3% 


- B 


2% 


- » 


1% 


-e 


0% 



«,'! 

(;'( 
>/« 

-0.02 1- is 

ft* 
4" 

"/I 

a 
K 

4 

6 



. 

6 



10 



15 
X 



20 



25 



30 



Figure 5.6: Diffusion barriers minus 
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layer as obtained in our model. Calcu- 
lations are shown for a single adatom 
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height. The misfit is chosen between 
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In order to facilitate the comparison of the results for different misfits we subtract 
E^ WL (e) from the activation barriers for diffusion jumps on the island resulting in AE a . 
Again we choose 24 particles as the base size of the exemplary island. Figure 5.6 shows 
results for a three layer high island. Two important features can be deduced: 

• For misfits e > the diffusion on top of the islands is in general slower than on 
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the wetting-layer (AE a > 0). The higher the misfit the slower becomes diffusion 
on the island compared to the diffusion on the the wetting-layer. 

• The diffusion bias towards the island center is more pronounced for higher misfits. 

In conclusion an increasing misfit should enhance island formation. As we will see a 
higher misfit indeed implies a larger island density and smaller mean island base sizes. 

5.4 Stranski-Krastanov-like growth scenario 

In our investigation of SK-like growth we follow now a scenario which is frequently 
studied in experiments [110,114]: in each simulation run a total number of 4ML ad- 
sorbate material is deposited at rates 0.5ML/s < Rd < 9.0ML/s. After the deposition 
is complete, a relaxation period with Rd = of about 10 7 steps follows, corresponding 
to a physical time on the order of 0.3s. Results have been obtained on average over at 
least 15 independent simulation runs for each data point. 

In our simulation we observe the complete scenario of SK-like growth as described 
above. Figures 5.7 and 5.8 show different stages of growth at T = 500K for R d = lML/s 
and Rd = 7ML/s, respectively. Both figures show a section of a L = 800 system for (top 
down) 1.5ML, 2.3ML, 3.0ML, 4.0ML coverage and after the mentioned relaxation time 
of 10 7 simulation steps. Ultimately, the formation of islands is driven by the relaxation 
of strain. Material within the 3d island and at its surface can adapt a lattice constant 
close to that of bulk adsorbate. On the contrary, particles in the wetting-layer are 
forced to adapt to smaller distances. Note that also the uppermost substrate layer is 
affected slightly by the adsorbate. 

During deposition monolayer islands located on the wetting-layer undergo a rapid 
transition to bilayer islands at a kinetically controlled critical wetting-layer thickness 
h* (first reported in [115]). We identify this transition as the 2d — 3d transition. The 
comparison of figures 5.7 and 5.8 already shows that h* increases with increasing de- 
position flux. For Rd = lML/s (see fig. 5.7) multilayer islands already start to form 
for a coverage of 1.5ML, whereas for Rd = 7ML/s (see fig. 5.8) the onset of island 
formation is delayed to about 2.3ML coverage. In particular from figure 5.8 it is obvi- 
ous that islands grow not only by capturing newly deposited particles, but also by the 
decomposition of the adsorbate layer. After the 2d — 3d transition the thickness of the 
wetting-layer decreases to a stationary value h c . Note from figure 5.8 that within our 
simulations it is possible for a very large island to split into well separated smaller ones. 

5.4.1 Kinetically controlled critical wetting-layer thickness 

For a systematic determination of h* we follow [114] and fit the density p of 3d islands 

as 

p = po (6 - Kf , (5.3) 
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Figure 5.7: Sections of a L = 800 system at T = 500fT and i? d = l.OML/s for (top 
down) 1.5ML, 2.3ML, 3.0ML, 4.0ML coverage and after 10 7 additional simulation 
steps. The darker the grey level of a particle, the bigger the average distance to its 
nearest neighbors of the same kind. The color bar gives a legend for the adsorbate 
particle distances in the downmost picture. 




Figure 5.8: Same as figure 5.7 but for = l.OML/s. Note that the onset of island 
formation is here delayed, resulting in h* ~ 2.3ML. 
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where 6 is the coverage with adsorbate material. The best fit for the data is obtained 
for po = 7 x 10~ 3 and a ~ 1.5. For example, in the case of InAs grown on GaAs one 
finds a ~ 1.76 [114]. Figure 5.9 shows the fit for T = 500K and different values of the 
deposition rate. 
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Figure 5.10 shows h* as a function of the deposition flux F for two temperatures. 
Increasing values of h* are observed with increasing flux and decreasing temperature. 
Qualitatively the same dependence of the kinetically controlled critical thickness on F 
and T was reported in [110] for InP/GaAs heteroepitaxial growth. 

This behavior of h* can be understood within our model: we consider fairly low 
diffusion barriers at high temperatures and suppress the ES effect explicitly. As a 
consequence, layer-by-layer growth is favored and second layer nucleation will play a 
minor role in the formation of mounds. The limiting effect which sets the characteristic 
time for the 2d — 3d transition is the upward diffusion of particles at the edge of existing 
islands. Clearly, the rate for this process will decrease drastically with decreasing 
temperature, because a high activation barrier E up has to be overcome. This rate is to 
be compared with the time scale set by the incoming flux. A high flux F will fill the 
layer before particles can perform an upward hop, and hence it will delay the 2d — 3d 
transition. 

These considerations, together with the results of [110], suggest a functional depen- 
dence of h* c on F and T of the form 



K = h 



R 



up 



(5.4) 



where R up is the diffusion rate computed from E up of upward hops close to the transition 
according to an Arrhenius law. Equation (5.4) is not expected to hold for very low fluxes 
or high temperatures where h* nearly coincides with the stationary value h c [110]. 

In our simulations we observe at the 2d — 3d transition a typical value of E up pa 
leV for upward hops. Using this value we have performed a non-linear fit according 
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Figure 5.10: Kinetically controlled 
critical thickness h* vs. the deposi- 
tion flux F for two temperatures. Both 
curves correspond to equation (5.4), 
with h Q and 7 obtained from the data 
for T = 500K, only. The error bars 
represent as the standard error of the 
simulation results. 



to equation (5.4) with the data for T = 500K. It yields very good agreement for 
h ~ 3ML,7 rs 0.2, cf. figure 5.10 (dashed lower curve). Setting the temperature to 
T = A80K in equation (5.4) with otherwise unchanged parameters yields the dashed- 
dotted upper curve in figure 5.10. The good agreement with the simulation data strongly 
supports our assumptions. Note that 7 is expected to depend strongly on material 
systems. Clearly our data does not allow for a precise determination of 7, however 
a positive value of 7 captures the essential F and T dependencies of h* like observed 
in [110]. 

5.4.2 Island properties 

In order to characterize the self-assembled islands that emerged during growth, we 
measure their base b given by the number of particles in the island bottom layer and 
their height h in monolayers. Figure 5.11 shows base size and height distribution for 
R d = l.OML/s and R d = 7.0ML/s, both at 500fT. For R d = l.OML/s the mean base 
size is 27 particles and the mean height is 10ML. For R d = 7.0ML/s the mean base 
size and mean height are 23 particles and 6 ML, respectively. 

As one can further conclude from figure 5.12 (left panel) the mean base size first 
decreases with increasing particle flux but becomes constant and independent of the 
temperature at higher values of the deposition rate. All results shown here are obtained 
at the end of the relaxation period with R d = 0. Whereas the mean values do not change 
significantly, fluctuations decrease with time in this phase. We also consider the density 
p of multilayer islands, i.e. their total number per particle in a substrate layer. Figure 
5.12 (right panel) shows p as a function of deposition rate and temperature. For low 
deposition rates the density rises with increasing R d while it decreases for increasing 
temperature. 

At higher values of the deposition rate the island density becomes constant. A simi- 
lar behavior of the island density as a function of R d and T was reported for InP/GaAs 
heteroepitaxial growth [110,126]. The constant, temperature independent region of the 
density is reminiscent of the behavior of the base size for high fluxes. The saturation 
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Figure 5.12: Left panel: Average base size < b > of multilayer islands vs. the 
deposition rate Rd for two temperatures. Right panel: Density of multilayer islands 
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standard error of the simulation results. 
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behavior further demonstrates the importance of upward hops vs. aggregation of de- 
posited particles on islands. The latter process would yield a continuous increase of the 
island density with the flux. 

5.4.3 Influence of the misfit 



Figure 5.13: Average base size < b > 
of multilayer islands vs. the misfit e at 
R d = A.5ML/s and T = 500K. 
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Since our above calculations indicate that the misfit has an important impact on 
the island formation, we finally take a look at the misfit dependence of the island size 
and density. From several simulation runs for values of the misfit 2% < e < 4% at 
Rd = 4.5ML/s and T = 500K we deduce that the base size at high deposition rates is 
mainly determined by the misfit and decreases with increasing misfit like b ps 0.9i£ -1 
(see fig. 5.13). Figure 5.14 shows corresponding sections from simulation runs for 
R d = A.bML/s and T = 500fT. In misfit heteroepitaxy due to the incompatibility of 
the substrate and the adsorbate lattice, e' 1 is the relevant length scale (cf. chapter 
4). The driving force for the misfit dependence of the island size is probably the slower 
diffusion on top of islands and the increased bias to the island center for higher values 
of e, cf. figure 5.6. We believe that the observed dependence of the base size on the 
misfit for high growth rates is reminiscent of this fact. High growth rates lead therefore 
to a narrow size distribution around a self-limiting, strain dependent island size. 

Note that in situations close to equilibrium one would expect an e-dependence 
b oc s~ 2 , see e.g. [15] for theoretical arguments. In principle, this regime should be 
attained within our model in the limit Rd — > where the observed size b and island 
height increases drastically, indeed. 

5.5 Conclusions and outlook 

In conclusion we present a model for the simulation of SK-like growth which is capable 
of reproducing various important phenomena observed in experimental studies. Our 
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Figure 5.14: Typical islands 

for (top down) e = 2%, e = 3% 
and e = 4% for deposition rate 
R d = A.bML/s and temper- 
ature T = 500fT. The pic- 
tures show sections of the same 
size from simulation runs with 
L = 800. 



work provides a fairly detailed and plausible picture of the SK growth mode. We have 
shown that a strong adsorbate-substrate interaction can be the cause for the formation 
of a stable wetting-layer due to the relatively slow diffusion of adatoms on the substrate. 

We have demonstrated that the appearance of a kinetically controlled critical wet- 
ting-layer thickness, like observed in various experimental works, can be explained 
within our model. We argued that high particle fluxes and low temperatures stabilize 
the metastable adsorbate layer. 

The formation of islands is traced back to the relaxation of strain in the adsorbate 
layer: this strain relaxation leads to a pronounced bias towards the island center on top 
of the multilayer island. In addition the diffusion is slower on top of the island than on 
the wetting-layer in our model. The island size is determined by the particle flux, the 
temperature and most importantly the misfit between adsorbate and substrate. For 
large enough particle fluxes, the island size becomes independent of the temperature 
and is controlled by the misfit only. 

Though we were able to gain a first insight into relevant mechanisms of self- 
assembled island formation various important questions still remain open. One inter- 
esting point is whether intermixing between substrate and adsorbate material stabilizes 
the wetting-layer and thus increases the critical wetting-layer thickness. To this end 
exchange diffusion in the bulk has to be considered in the model. Like for all other 
types of concerted moves it is not yet clear if this is conceptually possible within our 
method. 

Ultimately the method should be extended to 2 + 1 dimensions, more realistic empir- 
ical potentials and realistic lattice structures. Due to the high diffusion rates and low 
particle fluxes, which are necessary for the simulation of Stranski-Krastanov growth 
this is currently beyond computational means. 
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Chapter 6 

Simulation of multi-component 
growth 

In this chapter we focus on the analysis of multi-component growth by means of Monte 
Carlo simulations. To this end we study the submonolayer growth of two types of ad- 
sorbate particles on a given substrate. Our examinations are mainly motivated by ex- 
perimental studies, where a variety of material systems have been found which, though 
immiscible in the bulk, form stable alloy layers if deposited as a thin film on certain 
substrate materials. 

Though this phenomenon is witnessed for various material systems (for example 
CoAg/Ru(0001) [37-40], CoAg/Mo(110) [41], FeAg/Mo(110) [41], AgCu/Ru(0001) [42], 
PdAu/Ru(0001) [43]) they all have one thing in common: the types of deposited ad- 
sorbate particles show a misfit of opposite sign with the substrate. 

Indeed it was shown [31] that in systems dominated by an atomic size mismatch 
surface confined alloying is a possible strain relaxation mechanism (also see chapter 
1). This is due to the fact that mixing between the different types of material reduces 
the strain energy in the surface layer. Within enhanced models [127, 128] - taking the 
competition between chemical binding and elastic energy into account - it was shown by 
means of equilibrium simulations, that this competition can result in a striped structure 
of the film. Depending on the concentration of the adsorbate materials and the chemical 
binding between them, alloying and the formation of misfit dislocations are competing 
relaxation processes [39,127]. 

These theoretical findings are also in agreement with various experimental studies, 
where under certain growth conditions both adsorbate types form islands of a regular 
stripe structure. Furthermore the alloying can lead here to a change in the morphology 
of the islands' shape. In case of CoAg/Ru(0001) for example both adsorbate types 
form islands of compact shape when deposited alone on the substrate. But their co- 
deposition yields islands of pronounced ramification riddled with alternating veins of 
approximately constant width (see fig. 6.1). 

Although the observed striped islands show a remarkable similarity with the equi- 
librium simulation results it is not yet clear if strain effects are indeed the driving force 
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Figure 6.1: 5000 x 5000 nm scan- 
ning tunneling microscopy image of a 
dendritic CoAg island on Ru(0001) (by 
courtesy of R.Q. Hwang [37]). 



for the observed phenomena. As proposed in [37] also a difference between the chemical 
binding energy of equal and unequal adsorbate particles would both favor the observed 
vein structure and the dendritic island shape. 

In the following we will present a model which is capable of simulating multi- 
component growth in 2 + 1 dimensions, incorporating both strain and binding energy 
effects. We will first discuss the results of some equilibrium simulations and compare 
them to KMC simulation results afterwards. We will try to distinguish whether the 
phenomena observed in KMC simulations are due to elastic or binding interactions. For 
that reason we also propose a lattice gas KMC model which is mainly governed by an 
edge diffusion barrier and where strain effects are completely neglected. 



6.1 Simulation model 

Since for conceptional reasons phenomena like the formation of alternating vein struc- 
tures or the ramified island growth can not be mapped to 1 + 1 dimensions we have to 
extend our simulation method to 2 + 1 dimensions. In order to keep the computational 
effort acceptable we choose the simple cubic (sc) symmetry for our simulations, which 
has two important advantages: first one has to take only four possible in-plane tran- 
sition sites per particle into account, in comparison with six transition sites in case of, 
e.g., the fee lattice. And more importantly, since the sc lattice is not close-packed the 
number of particles to be considered during the relaxation processes is rather small. 

But the choice of the sc symmetry implies also an important problem: all of the 
experimental results presented in [37-43] are given for metals grown on a substrate of 
fec/hep symmetry. However, this should primarily affect the (compact) island shape: 
in case of the discussed materials one finds triangular shaped islands, whereas in the sc 
symmetry square shaped islands should be formed. 

In order to stabilize the sc lattice we adapt the method proposed in [11] and choose 

M" + '®-aeH)GH))*- ™ 
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as interaction potential between two particles i and j (see appendix A. 3 for details). 
In the following two kinds of pair-potentials are used during our simulations: The 
Lennard- Jones (LJ) potential discussed in appendix A.l 




12 



) 
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(6.2) 



and the Morse potential (see appendix A. 2) 




)-2). 



(6.3) 



In order to save computer time during energy calculations Uij is cut off for particle 
distances greater than r cut = 2ro, whereas for the calculation of the activation barriers 
the cut-off distance is set to 3tq. Both simplifications are perfectly justified since the 
used potentials decline fast towards zero for increasing particle distance. For both types 
Eij gives the depth of the potential for the interaction of two particles i and j at the 
equilibrium distance r with r oc a^. In the following the interaction strength between 
two particles of the substrate is given by = E s and = as = 1. E^ = Ea, <Jij = a a 
and E^ = Eb, <Jij = ob are chosen for the A-A and B-B interaction, respectively. Here 
A and B denotes the two adsorbate species and S labels the substrate material. For 
the interaction of adsorbate particles of type X = A, B with the substrate we set 
= E xs = \/E x E s and a tj = a xs = l/2(o* + <r s )- Likewise, = E AB = \/E A E B 
and Uij = oab = 1/2(o"a + a B ) holds for the interaction between A and B adsorbate 
particles. The misfit e (e > 0) is applied in a symmetric way to the system: 



The additional parameter a in equation (6.3) determines the steepness of the Morse 
potential around its minimum. It is used with three different values of the parameter 
a: a = 5.0, 5.5 and 6.0. 

The potential depths E^ are chosen in such a way that they meet two demands: on 
the one hand the ratio between E s and E A) E B is kept fixed for all potentials: 



On the other hand the barrier for diffusion on plain substrate E a su b in the case of 
homoepitaxy (e = 0) should become roughly the same value for all used potentials to 
facilitate the comparison of the results. We choose Eg here in such a way that in the 
case of homoepitaxy (e = 0) E a:Sub ps 0.37eV - a typical value for self-diffusion barriers 
of metals (see e.g. [62,63, 129, 130]). The resulting E s for the different potentials are 
shown in table 6.1. 



a A = l-e 
a B = l + e. 



(6.4) 
(6.5) 



Ea — E B — —E$- 



(6.6) 
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potential 


LJ 


a=5.0 


a=5.5 


a=6.0 


E s /eV 


3.0 


3.0 


2.814 


2.70 



Table 6.1: The substrate-substrate interaction strength E s for the used pair- 
potentials. The substrate-adsorbate and adsorbate-adsorbate interactions are then 
given according to equations (6.6) and (6.9), respectively. 

6.2 Equilibrium simulations 

In order to derive the influence of the potential depth for the A-B interaction we first 
carry out canonical equilibrium simulations at full coverage of the substrate, where the 
number of A and B particles is kept fixed. The substrate is given as an 100 x 100 
crystal of 6 layers height. Periodic boundary conditions are applied in the x- and y- 
direction. At the beginning of each simulation run the substrate is randomly covered 
with adsorbate particles at a certain ratio ua/ub between the number of A and B 
particles. Then the system approaches thermal equilibrium at temperature T by means 
of a rejection-free method (see [4]). 

6.2.1 Rejection-free method 

Since the simulations are carried out for a range of the misfit where even at full coverage 
of the substrate misfit dislocation are not observed each particle can be allocated at 
a certain site of the 100 x 100 square lattice. In order to realize a fixed number of 
each particle species like required in canonical simulations we follow here a method for 
the rejection-free canonical equilibrium simulation proposed by Ahr [4,131]. In each 
event an A particle exchanges its binding site with a B particle. We choose a nonlocal 
dynamics where the range of particle jumps is unlimited. This yields considerably faster 
equilibration compared to local dynamics. For simplicity, we permit only exchange 
jumps of particles which are more than r cut away from each other, i.e. exchange of two 
particles within this range is forbidden. Such processes would complicate the calculation 
of the configuration energies. 

If now an A particle located on site % of the square lattice exchanges its binding site 
with a B particle of site j the energy difference between the final and the initial state 
is given by AH = AHj — AHi, where AH X = H X (A) — H X (B) is the energy difference 
of the system with site x occupied with an A and B particle. H X (A) and H X (B) are 
calculated in a local way: an A particle is set to the site x and the particles of the system 
within a radius r cut around this site are allowed to relax locally. The local energy is 
registered as H X (A). Conversely, the same procedure is performed for an B particle at 
site x in order to calculate H X (B). The rates 

Af/j-A Hj 

ri_j = e 2fer (6.7) 
fulfill the detailed balance condition. Then the probability for an A particle on site i 
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to exchange its site with a B particle at site j factorizes, i.e. 

^ = P^ B pt A where e.g. p^ B = (6.8) 

Here p A ^ B , f' A ~^ B give probability and rate for a change from A to B particle at site 
x, respectively. Analogous p B ~* A , f' B ^ A are given for a change from B to A particle at 
site x. The rate r A ^ B is given by r A ^ B = exp (AH x /2kT) if site x is occupied by an 
A particle and zero otherwise. Conversely, r B ^ A = exp (—AH x /2kT) on sites occupied 
by B particles and zero on sites occupied by A particles. 

Due to this factorization property one can proceed now in two steps: In the first 
step, using a binary search tree a site % occupied with an A particle is drawn with the 
probability p A ~^ B ■ Then a B occupied site j is selected with probability p B ^ A , using a 
second search tree. If the distance between both sites is greater than r cut the exchange 
jump is performed and the rates of all affected events are updated. Otherwise, the event 
is rejected and the system remains unchanged. Since the number of rejected events is 
small for large systems, the loss of speed can be neglected. 

In order to avoid artificial strain accumulation due to the local relaxation for the 
calculation of AH X the system is globally relaxed after a fixed number of simulation 
steps (here 5000), all rates are newly calculated and the search trees are updated ac- 
cordingly. The system's total energy per particle E tot is registered after each global 
relaxation. All simulation runs are halted after 20 global relaxation events (i.e. after 
10 5 simulation steps). 



6.2.2 Influence of misfit and binding energy 




Figure 6.2: 80 x 80 sections for the Lennard- Jones potential with Eab = at T = 
250K. The misfits is (from left to the right) e = 0%, e = 3%, e = 4.5% and e = 5.5%. 
The bigger B particles are shown in light gray. 

In order to detect the influence of misfit and binding energy we first present some 
results for the cubic Lennard- Jones potential (eq. (6.1), (6.2)). Figure 6.2 shows typical 
sections for simulation runs with Eab = at T = 250/T. Since there is no interaction 
between A and B particles for the homoepitaxial case (e = 0%) a fast demixing of the 
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two particle types is observed. A straight line in the x- or ^/-direction separates both 
areas minimizing the energy loss of the system. 

From the sample for e = 3% it appears, that a straight boundary between the 
both particle areas is less favorable. This becomes more obvious for higher misfits (e.g. 
e = 4.5%, e = 5.5%) where several clusters of both particles types are formed and 
the boundaries between A and B domains are preferentially aligned in the < 11 > 
directions. 

The reason for the observed e-dependence is quite clear: the higher the misfit the 
less favorable become extended areas of the same particle type. The strain energy in 
such big clusters even exceeds the energy lost due to the enlarged domain wall between 
regions of A und B particles. However, even though the misfit causes separation between 
the particle types at Eab = 0, no regularity or ordering of the structures (like the 
alternating veins mentioned before) is observed. 

The situation changes completely for the case of Eab > (see fig. 6.3). Now 




Figure 6.3: Snapshots for the Lennard-Jones potential at T = 250K with Eab = 
0.6E A , E AB = 0-8E A , E AB = 0-9E A , E AB = 1.0E A (from the left to the right) and 
e = 4.5%, e = 5.0%, e = 5.5% (top down). The panels for E AB = 1.0E A show 40 x 40 
sections of the system, the remaining panels show 80 x 80 sections. 



a regular arrangement of alternating A and B stripes, preferentially in the < 11 > 
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directions is visible. As already known from other atomistic models with size mismatch 
[31, 128] the competition between binding of the particles and strain energy is the 
cause for these regular patterns. Furthermore with, increasing Eab and increasing 
misfit the stripes become thinner and more regular in size and shape. For the case 
Eab = Ea = E b the system approaches a checkered state, i.e. a stripe width of one. 

The alignment of the stripes in the < 11 > direction - already visible for the E A b = 
case (see fig. 6.2) - is here due to the cubic symmetry of the used potential: both particle 
types try to reach their preferred stripe width / in each lattice direction (x and y). Note, 
that the used cubic form of the potential (eq. (6.1)) has only a weak interaction in the 
< 11 > direction (also see fig. A. 3 in the appendix A). 
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Figure 6.4: (a) Stripe width for Eab = 0.6Eab and particle concentrations t)a = 
T] B = 0.5 as a function of the misfit, (b) Stripe width for E A b = 0.9-EU.b, £ = 5% as a 
function of the A particle concentration r\ A {j]b = 1 — Va, consequently). Temperature 
is T = 250i^ for both pictures. Each value is obtained by averaging over 3 independent 
simulation runs. The errorbars are given by the standard deviation. 



Figure 6.4(a) shows the width I of A and B stripes for Eab = 0.6Ea as a function 
of the misfit. Since the concentrations tja — r)B — 0.5 of A, B particles are the same 
A and B stripes have about the same width, whereas the situation changes completely 
when i]a 7^ t]b- As figure 6.4(b) shows for Eab — 0-9Ea and e = 5% the stripe width 
increases with increasing concentration of the particle type. It is noticeable that the 
bigger B particles form thinner stripes at high B concentration than the smaller A 
particles at high A concentration. This is due to the asymmetric pair-potential, which 
is steeper in compression than in tension and thus (compressed) B stripes are a little 
more restricted in their width than A stripes. 



6.2.3 Temperature dependence 

We take now a short look at the temperature dependence of the stripe structures. To 
this end simulations for the Lennard- Jones potential are performed for T = 500i^. The 
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simulations were halted after the same number of simulation steps as for T = 250K. A 
comparison of samples at T = 250K and T = 500K shows, that in regions of smaller 
binding energies E AB < 0.6E A the higher temperature yields widened stripes. This 



Figure 6.5: Snapshots (80 x 80 sections) 
for the Lennard- Jones potential at T = 
500fT at e = 4.5% for E AB = 0.6E A (left 
panel) and E AB = 0.9E A (right panel). 



is due to the fact that the system approaches a deeper value of E tot for the higher 
temperature in shorter times (also see fig. 6.6), i.e. a smaller number of simulation 
steps is required as for T = 250K. For higher values of E AB the high temperature only 
results in rather noisy stripes, but of about the same width as for T = 250K. 



Figure 6.6: Total energy of the sys- 
tem per particle E tot vs. the number 
of simulation steps for T = 250K and 
T = 500JT. E AB is set to O.QE A and 
the misfit is e = 4.5%. Each data point 
is obtained by averaging over 3 inde- 
pendent simulation runs. The errorbars 
represent the standard deviation. 

■ .. . 10*5 

simulation steps 



6.2.4 Influence of the interaction potential 

With otherwise unchanged parameters the simulations are now performed for the a = 
6.0 Morse potential, which is steeper in both - compression and tension - than the 
Lennard- Jones potential used before (also see appendices A.1,A.2). The simulation 
temperature in set to T = 250K. However, as figure 6.7 shows Lennard- Jones and 
Morse a = 6.0 potential yield quite similar results: again the competition between strain 
and binding energy causes alternating stripes of decreasing width with increasing e. Due 
to the cubic symmetry the stripes are again solely aligned in the < 11 > direction. 

As figure 6.8 points out the main difference one observes is that for the same mis- 
fit and E AB < §.&E A the stripes for the Morse potential are systematically thicker. 
Whereas at higher values of E AB the mean stripe width is nearly identical for both 
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Figure 6.7: Snapshots for the Morse a = 6.0 potential at T = 250fT with E AB = 0.6E A , 
E AB = 0.8E A , E AB = 0.9E A , E AB = 1.0E A (from the left to the right) at e = 5.5%. 
The panel for E AB = 1.0E A shows a 40 x 40 sections of the system, the remaining 
panels show 80 x 80 sections. 
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Figure 6.8: Width / of B stripes 
as a function of E AB for different 
values of e. Each data point is ob- 
tained by averaging over 3 indepen- 
dent simulation runs. The error- 
bars are given by the standard de- 
viation. 
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potentials at a given misfit. However, even at values E AB < 0.QE A the deviation are 
small compared to the influence of the particle concentration or the temperature on the 
stripe width. 

In conclusion we learn from the equilibrium simulations, that for the heteroepitaxial 
case (e > 0) the state of minimum free energy is not given by a complete separation 
of the different particle types, but rather splits the system in several clusters of both 
species. A non-vanishing binding energy E AB > between both particles types together 
with e > yields regular patterns of alternating stripes. The width of these stripes is 
above all controlled by the value of e together with the binding energy. The competition 
between strain and binding energy results in structures not unlike the ones observed 
experimentally. However, the question is now if also growth simulations of a system, 
which is governed by these competition, yield a similar morphology. 
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6.3 KMC simulations 



In the following we discuss results obtained in off-lattice KMC simulations. Just like in 
the equilibrium simulations growth takes place on a 100 x 100 substrate of 6 layers height 
with periodic boundary conditions applied in the x- and ^-direction. For all simulation 
runs the deposition rate for both types of particles is set to 0.005ML/s resulting in an 
overall deposition rate of = O.OlML/s. The simulations are halted when half of 
the substrate is covered with adsorbate particles. Since we are only interested in the 
submonolayer regime here, we disregard second layer nucleation, i.e. particles which 
are deposited onto other particles will be ignored. This implies the existence of a 
rather high Ehrlich-Schwoebel barrier for diffusion across step edges. Note that for the 
same reason jumps of particles onto others are suppressed. Unless otherwise mentioned 
the temperature is set to T = 500i^ in the following. The same attempt frequency 
uo = 10 12 Hz is considered for all diffusion processes. 



Figure 6.9: Schematic representa- 
tion of the criterion for barrier cal- 
culation. Only if a adsorbate parti- 
cle feels other particles within 
dius of 3tq barriers are calculated. 
Otherwise a stored value for diffu- 
sion on plain substrate is selected. 



For the interaction between A and B particles we choose 

E AB = 0.6E A = 0.6E B = 0.10E S . (6.9) 

- a value well between the ones investigated during the equilibrium simulations, where 
rather thick stripes are formed and the influence of the misfit should be clearly observ- 
able. Furthermore, on the basis of the equilibrium simulation results, we expect to find 
a marked dependence on the choice of the potential for this interaction strength. Here 
again Ea = E B are given according to equation (6.6) and Es is chosen from table 6.1 
for the different potentials. 

Note that this choice of the potential depth yields a higher barrier for edge diffusion 
than for diffusion on plain substrate in our simulations. However, the barrier for edge 
diffusion is still smaller than that for detachment from the edge. So particles attached 
to an island edge are more likely to diffuse there than to detach. This is of particular 
importance since we focus here on phenomena, where edge diffusion is supposed to 
have a strong impact. Note also that for the cubic lattice (eq. (6.1)) diagonal diffusion 
jumps can be neglected since as figure 1.2 shows the barrier for a diagonal jump is here 
represented by a maximum in the PES. 
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The realization of the simulation algorithm is as described in chapter 3. The simu- 
lations are carried out for a range of misfits in which even at full coverage dislocations 
are not observed and each particle is allocated at a certain site of the square lattice. 
We can therefore take advantage of the lattice based method proposed in chapter 3. 

Since we simulate submonolayer growth we are able to adopt a further simplification 
to the method. Due to the cut-off distance for barrier calculation, at a given misfit 
the diffusion barrier for an adsorbate particle on plain substrate does not change as 
long as no other adsorbate particle is within a radius of 3r (see also fig. 6.9). For this 
case stored values can be used for the diffusion barrier of A and B particles on plain 
substrate, introducing only very little error. Additionally to the barrier calculation also 
the local relaxation can be omitted here. Especially in the early stages of growth this 
accelerates the computations a lot. Preceding simulation runs with full calculation of all 
diffusion barriers showed that these fixed barriers for diffusion steps on plain substrate 
have no impact on the results of our simulations. 
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Figure 6.10: Snapshots for the Morse a = 5.0, a = 5.5 and 6.0 potential (top down) 
for different values of e. The bigger B particles are shown in light gray. 
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Figure 6.11: (a) Ratio A between perimeter particles and total number of particles 
in the big B clusters for the used potentials, (b) The number of perimeter particles 
divided by the square root of deposited particles T vs. e. Each value is obtained by 
averaging over 10 independent simulation runs. The errorbars are given by the standard 
deviation. 



6.3.1 Influence of misfit and potential 

We present now results on the influence of the value of e and the used potential at a tem- 
perature T = 500i^. Preliminary simulation runs showed that under the same growth 
conditions both particle types form rectangular shaped islands if they are deposited 
alone on the substrate. The situation changes completely in the case of co-deposition: 
Figure 6.10 shows snapshots of simulation runs for the Morse a = 5.0, a = 5.5 and 
a = 6.0 potential for 4 different values of e. 

These structures are exemplary for all simulation results: the B particles (shown in 
light gray) assemble into a few big clusters. With increasing misfit the branches of these 
clusters become thinner and of more uniform width. The A particles surround these 
branches without showing a similar shape. It is also seen from figure 6.10 that with 
increasing misfit the ramification of the structure as a whole increases. This is clearly 
related to the restricted width of the B stripes: a B particle rather attaches to the thin 
end of a stripe. That means the thinner the stripes the faster is the outwards growth 
of the light gray branches, leading to an increasing ramification of the structure. 

At a given misfit the B branches are the thinner the smaller the value of a in the 
Morse potential is. Consequently, at a given misfit the island-ramification is more 
pronounced for a = 5.0 than for a = 6.0. This is in agreement with the equilibrium 
simulations where a steeper potential yields thicker stripes. 

In order to quantify the made observations we calculate for each connected cluster 
of B particles the ratio A between its perimeter length and its volume. This is done by 
counting the number of perimeter particles together with the total number of particles 
in the same cluster. We take only the backbone of the structures into account and 
neglect smaller clusters (< 700 particles). The ratio A should then give a measure 
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for the average thickness of this cluster (see fig. 6.11(a)). For example, for a rather 
thin cluster most of its particles sit at the edge and therefore A should be close to 1, 
whereas with increasing cluster thickness A should decrease towards 0. In addition we 
measure T given by the number of particles with less then 4 nearest neighbors divided 
by the square root of the number of deposited particles (see fig. 6.11(b)). T provides a 
measure for the length of the structure's perimeter and therefore the ramification. For 
one perfect quadratic island on the substrate it results to T 4, whereas an increasing 
T indicates roughening of the island's shape. The correlation between A and T is clearly 
observable for all used potentials: A increases with increasing misfit indicating thinner 
B clusters. Simultaneously the ramification increases. The formation of B branches of 
well-defined thickness is a common phenomenon for the used pair-potentials. 

6.3.2 Influence of the temperature 






Figure 6.12: Snapshots for the Morse a = 6.0 potential, e = 6.5% at T = 400/T, 
T = A50K, T = 525K and T = 550fT (from the left to the right). 



One could now argue that the observed ramification of the islands is due to tem- 
perature effects: i.e. the used temperature is high enough for the formation of cubic 
clusters if a single particle type is deposited, but enlarged edge diffusion barriers in the 
co-deposition regime cause dendritic growth. 

In order to investigate the temperature dependence of island growth we simulate for 
the Lennard- Jones at e — 5.0% and the a = 6.0 Morse potential at e — 6.5% growth 
for different temperatures. As shown above for the given parameters strongly ramified 
islands grow at T = 500K. Figure 6.12 shows some snapshots for 400-K" < T < 550-K" 
for the Morse potential. Similar results are obtained for the Lennard- Jones potential. 
At T = 400-fT (fig. 6.12, left most panel) two islands emerge due to the reduced 
diffusion length. Both show frayed edges and rather thin and disordered B stripes. 
With increasing temperature the B stripes become wider and more regular in shape, 
the island edges become smoother. Also A and T as functions of T (see fig. 6. 13(a), (b)) 
exhibit this behavior. It is remarkable here, that the ramification T does not decreases 
monotonously with increasing temperature (as one would expect), but passes for both 
potentials through a minimum at T 475K. Then T slowly increases again with 
increasing T. This confirms that the observed ramification is not due to the low growth 
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T/K T/K 

Figure 6.13: A (a) and T (b) as functions of the temperature for the Lennard- Jones 
and Morse a = 6.0 potential at e = 5.0%, e = 6.5% respectively. 

temperature. The enhanced mobility of the particles causes a more distinct separation 
of the two particle types, resulting in more regular B stripes. As figure 6.13(a) shows 
the width of the B stripes approaches a constant value for the high temperature region. 
Further note, that in the region of high temperatures - like in the equilibrium simulations 
- nearly all B clusters are aligned in the < 11 > directions in order to achieve the 
energetically most favorable arrangement of particles. 

6.4 Lattice gas simulations 

The question now is in which way the observed branches are related to the stripe 
structures found in the equilibrium simulations. Figure 6.14(a) shows the activation 
barrier for a B particle diffusing near an A-B interface for e = 4% and E^b — 0.6EU. 
The weaker A-B interaction causes an extra step edge diffusion barrier for the jump 
from the B to the A region. As already mentioned this diffusion barrier is believed to 
favor the formation of alternating stripes. We discuss now by means of a latticed based 
method how such a diffusion barrier influences the multi-component growth. Of special 
interest is here if the stripe formation and the island morphology, as observed in our 
off-lattice simulations, can be explained within such a model. 

6.4.1 Simulation Model 

We consider now a lattice model with the two adsorbate species growing on a different 
planar substrate which provides a square lattice with 150 x 150 adsorption sites (x,y). 
Unlike in the off-lattice simulation where a particle interacts with all particles within 
the range of the potential, A and B particles interact now only with their lateral nearest 
neighbors through attractive two-particle interactions E AA , E BB and E AB . The dif- 
fusion of adatoms on the surface is described by thermally activated nearest neighbor 
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Figure 6.14: (a) Barriers for diffusion of a B particle (light gray) from the right to the 
left, near an A-B interface. The values are given for the Lennard-Jones potential at a 
misfit e = 4%. Note the increased barrier for a diffusion jump towards the A region, 
(b) Diffusion barriers for the diffusion process of a B particle in the lattice gas model. 
Note that the particle has to overcome an additional barrier E° — E AB > for a jump 
to the A region. 



hopping processes with Arrhenius rates Ri = z/oexp(— E a> i/kT), where we again choose 
uq = 10 12 Hz as value for the attempt frequency. Unless otherwise indicated the tem- 
perature T is again set to For the activation energy E a ^ of the diffusion event i 
we use Kawasaki type energy barriers [80]: 

Eaj, = max{B x , B x + AE X }. (6.10) 

Here B gives the diffusion barrier for a free particle of type X = A or B on the 
substrate. AE X denotes the change in the total energy caused by the event which is 
determined by a simple bond-counting scheme from the NN particle-interactions. For 
example, for a diffusion event of an A particle we obtain 

AE A = An AA E AA + An AB E AB , (6.11) 

where An and An AB count the difference between the number of A-A and A-B bonds 
before and after the diffusion step. To keep the number of parameters manageable we 
assume equal diffusion barriers for free A and B particles, B A = B B = B°, and also the 
strength of A-A and B-B bonds shall be the same: E AA = E BB = E°. The interaction 
between particles of different types is assumed to be weaker than between those of the 
same type: E AB < E°. 

Figure 6.14(b) shows activation energies for exemplary diffusion processes of a free 
B particle and B particles at the boundary of a small A-B cluster. For A particles 
the picture is essentially the same, except that the activation energies for crossing the 
A-B interface have to be exchanged, as well as those for detachment from A and B 
step edges. Since E AB < E° one reads from figure 6.14(b) that the difference between 
E AB and E° has two main consequences. First, an A particle diffusing along a step 
edge made up of A particles faces an enhanced diffusion barrier B° + E° — E AB > B° 
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when it attempts to cross an A-B interface. The same happens to a B particle coming 
from the other side. Thus, A and B particles diffusing along step edges are likely to be 
reflected at A-B interfaces. Second, the barrier for detachment of an A particle from a 
step edge made up of B particles is lower than that for detachment from an A step edge 
and vice versa. This reproduces basically the influence of a weaker A-B interaction in 
the off-lattice simulations, disregarding though all influences of strain and long range 
interactions. 

6.4.2 Influence of the step edge barrier 

We first investigate the influence of the binding energy E AB between A and B particles 
on the morphology of growing films. Therefore we fix B° = 0.37eV and E° = 0.51eV 
- this reproduces roughly the homoepitaxy (e = 0) barriers for diffusion on planar 
substrate and detachment from an island edge as measured in the off-lattice simulations. 
E AB is varied between 0.31£° and 0.71£°. 

Island geometry 

Following the off-lattice simulation for all simulation runs the deposition rate for both 
types of particles is set to 0.005ML/s resulting in a overall deposition rate of Rd = 
O.OlML/s. When the total adsorbate coverage has reached 0.5ML the particle fluxes 
are switched off and the simulation is halted. 

Figure 6.15 shows snapshots of example configurations obtained after the end of 
simulation runs for four different values of the binding energy E AB . 




Figure 6.15: Islands obtained by simultaneous deposition of A and B particles for 
different binding energies between A and B particles (from left to the right): E AB = 
0.71£°, E AB = 0.59£°, E AB = 0.51£° and E AB = 0.39£°. The system size is 150 x 150 
and the total coverage is 9 = 0.5ML. The figures are by courtesy of Th. Volkmann. 



For all values of E one observes compact island shapes with the island boundaries 
roughly parallel to the lattice directions. The weaker binding energy between A and 
B particles leads to an aggregation of particles of the same type in clusters which can 
be characterized as stripes. While for the higher value of E AB these stripes are rather 
thin and show a considerable degree of irregular intermixing for lower values of E AB 
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the stripes are both much thicker and there is a tendency for them to stretch outwards. 
One also sees that at a certain stage of the island growth a stripe of one particle type 
may become wide enough for particles of the other type to form a stable nucleus within 
this stripe, thus leading to a branch-like structure. The occurrence of the stripe-like 
structures is here a pure kinetic effect. 

Step geometry 

Due to the island topology chosen we have a fourfold symmetry of the stripe-structures. 
In the following we focus on the stripe formation in only one growth direction. Therefore 
we choose a different growth topology for our simulations where the growth proceeds 
from a step edge. We assume here that our system represents one particular terrace of 
a vicinal surface with elongated terraces of constant width. Therefore we use periodic 
boundary conditions in the direction parallel to the step edges (x-direction). In the 
perpendicular y-direction the system is bounded by the lower part of a step edge and 
the upper part of the next step edge. While adsorbate particles of both types may 
attach to the lower part of the step edge with the same binding energy E° as to next 
neighbors of the same type, they become reflected at the upper part of the step edge, 
again representing an infinite Ehrlich-Schwoebel barrier for interlayer jumps. 

We mostly studied systems with 256 x 100 and 512 x 100 lattice sites. The combina- 
tion of the specific boundary conditions in the y-direction together with the values used 
for the activation energies guarantees that for the given growth temperature T = 500i^ 
and system sizes adsorbate particles deposited on the substrate will reach the step edge 
at y — 0. Thus terrace nucleation is suppressed and growth proceeds from the step 
edge. All other simulation parameters remain unchanged. 

Figure 6.16 (left panel) shows snapshots of system configurations after the end of 
simulation runs for various values of E AB . The difference between E AB and E° leads 
to a separation of A and B particles in stripe-like clusters. Figure 6.16 (left panel) 
shows that for a decreasing value of E AB the average thickness of the stripes tends to 
increase. While for a low value of E AB A and B cluster are well separated they become 
more and more intertwined with increasing E AB . Both observations can be explained 
by the fact that with decreasing E AB it is more favorable for A and B particles to attach 
to particles of the same type. We notice also that there is no clear orientation of the 
stripes, for example parallel to the y-direction. In contrast stripes of one particle type 
may very well grow sideways as can best be seen in the downmost sample. 

We now have a closer look at the influence of E AB on the stripe width. Therefore 
we determine for each connected cluster of A particles the ratio A between its perimeter 
length and its volume. For each value of E AB we did 10 independent simulation runs 
and averaged over the occurring A clusters. To get a better statistics all clusters were 
initially sorted by their size and very small clusters (< 0.2 x the mean cluster size), for 
which A 1 holds, were omitted from the averaging. 

Figure 6.16 (right panel) shows the dependence of A on the binding energy E AB . As 
one sees A increases monotonously with increasing E AB confirming that the stripe width 



85 







J '^^^^^^^^^^^ ^^^^^^^^^^^^ ^^^^ ^^^^^^^^^ 

















0.8 
0.7 
0.6 
0.5 
<0.4 
0.3 
0.2 
0.1 



0.3 



0.4 



0.5 

E AB /E° 



0.6 



0.7 



0.8 



Figure 6.16: Left panel: simulation results for E AB /E° = 0.63,0.55,0.47,0.39,0.31 
(from top to bottom). A particles appear in dark gray, B particles in light gray. Right 
panel: dependence of A on E AB /E°. Each value is obtained by averaging over the A 
clusters of 10 independent simulation runs. The errorbars are given by the standard 
deviation. The shown results are by courtesy of Th. Volkmann. 



decreases with the difference in binding energies becoming smaller. For E AB /E° = 0.31 
we obtain A pa 0.08 which is already close to the value which we would expect for the 
limit E AB /E° — > when only one A cluster with straight edges was present. When 
E AB /E° approaches 1, A also should tend to 1 because A and B particles should then 
be perfectly mixed, since there is no longer a preference for a particle to stick to its own 
kind. This means that most of the cluster contain (9(1) particles which then results in 
Awl. 

We conclude from our lattice gas simulations, that the step edge barrier indeed 
gives reason for stripe formation. Similar to the equilibrium simulations the width 
of the stripes can be controlled by adjusting the binding energy between A and B 
particles. However as figure 6.15 shows neither asymmetries between A and B clusters 
nor a ramification of the islands is observed here. Certainly this is not surprising since 
in the lattice gas simulations A and B particles are treated in a symmetric way, whereas 
in the off-lattice simulations the different sign of the misfits causes different diffusion 
barriers for A and B particles, respectively. 



6.4.3 Enhanced lattice gas model 

For that reason we propose now an enhanced lattice model taking basic differences of 
both particle types into account. This is done in order to determine whether a simple 
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misfit dependence of the diffusion barriers could lead to the observed results within a 
lattice model. For example in the off-lattice case for e > the substrate diffusion for 
the bigger B particles is always faster than for the smaller A particles (cf. chapter 1). 

Furthermore in the off-lattice method the barriers for edge diffusion are higher 
than the substrate diffusion barriers. This could also give rise to a ramified island 
morphology. Therefore we extract the barriers for free diffusion on the substrate as 
well as averaged values for edge diffusion and detachment for a fixed island size (see 
also fig. 6.17) as a function of the misfit. Theses barriers are then used in the lattice 
model as parameters. The thus modified lattice model incorporates the basic misfit 
dependence of the diffusion barriers. But effects of the long-range interaction, like the 
reduced barriers for jumps towards an island still have to be neglected here. 
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Figure 6.17: Barriers for in-plane diffu- 
sion of a B particle (light gray) from the 
right to the left, near an A-B interface. 
The values are given for the Lennard- 
Jones potential at a misfit e = 4%. Note 
that for the off-lattice barriers the parti- 
cle feels the influence of the island due to 
the range of the potential before attaching 
to it. This results in a rather low barrier 
for the jump towards the island. 



Figure 6.18 shows a comparison between the lattice model and the off-lattice sim- 
ulation for the Lennard-Jones potential. Simular results are obtained by fitting the 
barriers for the Morse potentials to the lattice model. As expected, the islands for both 






Figure 6.18: Snapshots for the Lennard-Jones potential for the enhanced lattice and 
the off-lattice model. The panels show (from left to the right) lattice, off-lattice results 
at e = and lattice, off-lattice results for e = 5%. The lattice gas results are by 
courtesy of Th. Volkmann. 



models look similar in the case of zero misfit. However for e = 5% the results for both 
models seem to have little in common. For the lattice model the separation of A and 
B regions is more pronounced as for the e = parameters but no size limitation of 
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the stripes and therefore no ramification is observable here. To confirm this we have 
measured the ramification V for the lattice model. Figure 6.19 shows this in comparison 
to the values obtained from the off-lattice simulations. Only at e = both curves col- 
lapse. From this examinations it is quite clear that the basic differences of the diffusion 
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Figure 6.19: Comparison of the island- 
ramification in case of lattice and off- 
lattice simulations. Each value is obtained 
by averaging over 10 independent simula- 
tion runs. In case of off-lattice simula- 
tions the errorbars are given by the stan- 
dard deviation. For the lattice simulation 
errorbars are smaller than the symbols. 



barriers of both particle types neither cause the width restriction of the B branches nor 
the ramification of the islands with growing misfit. 



Figure 6.20: Snapshots for simulation 
runs with (left panel) and without (right 
panel) the reduced barrier for jumps to- 
ward island. Both samples are given for 
the a = 6.0 Morse potential at e = 6.5%. 



From additional off-lattice simulations, where the reduced barrier for a jump towards 
an island (cf. fig. 6.17) is suppressed we find that - though the resulting islands are less 
ramified - the width of the B branches stays the same (see fig. 6.20). The reduction 
of the ramification is due to a higher mobility of the particles: once a particle detaches 
from an island it has the same probability for jumps towards the island as away from 
it. The capturing of diffusing particles by islands is less pronounced and therefore the 
particles are more uniformly distributed around the island. 

In conclusion we learn from the enhanced lattice model with fitted diffusion barriers 
that neither the ramification nor the width restriction of the B stripes is due to the basic 
differences of the diffusion barriers of the both particle species. However, there has to 
be a kinetic reason for the obvious difference between both species. For that reason 
we take in the following an closer look on the diffusion behavior of particle attached to 
island edges. 
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Figure 6.21: Barriers for diffusion of an A particle near and on an A edge (a) and for 
a B particle diffusing near and on a B edge (b). The barriers were calculated for the 
Morse a = 6.0 potential, stripe width Z = 11 and misfit e = 5%. 



6.5 Diffusion bias for the edge diffusion 

By taking in our lattice model only averaged diffusion barriers into account we neglected 
the diffusion bias for edge diffusion. This bias is thereby due to the same reason as 
explained in chapters 1 and 5 for the diffusion on an island: near the center of an edge 
the particles are bound with the substrate lattice constant, towards the rims of the 
island they can approach their own lattice constant. Figure 6.21 shows the resulting 
activation energy E a for a particle diffusing near and on an A (e < 0) or B (e > 0) 
island edge. Here, the diffusing particle is of the same type as the island. However, the 
direction of the diffusion bias only depends on the type of the island, whereas the value 
of the barriers changes according to the interaction strength Eab between the particle 
types. 
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Figure 6.22: Relative difference 
AE e dge between the minimum and the 
maximum value of the edge diffusion 
barrier (width / = 11) as a function of 
the misfit e. 
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One reads from figure 6.21(a) that in case of an A edge the particle diffusion is biased 
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towards the rims of the islands, whereas for B edges particles are more likely to jump 
towards the center of the edge. This explains qualitatively the different morphology 
of A and B clusters in our simulations. The diffusion bias towards the center of a B 
edge (see fig. 6.21(b)) assists the formation of width restricted B stripes. The increased 
probability for particles to be found near the rims of A edges causes the observed 
broader structures. 

Within this framework also the dependency of the B stripes width on e and the 
used potential is understood. Figure 6.22 shows the relative difference AE e d ge between 
the minimum and the maximum value of the barrier for diffusion in one direction of a 
B particle attached to a B edge. The values are taken for a fixed stripe width (I = 11). 
For the considered potentials increases AE edge with increasing misfit, i.e. the influence 
of the bias becomes stronger. It is also seen from figure 6.22 that at a given misfit 
the bias is the stronger the shallower the potential is. Both observations correspond 
well with the behavior of the stripe width as a function of misfit and potential (see fig. 
6.11). 

6.6 Conclusions and outlook 

In this chapter we have demonstrated that our off-lattice method is capable of simulat- 
ing multi-component growth in 2 + 1 dimensions. It was shown that surface confined 
alloying of the two adsorbate species is indeed a possible strain relaxation mechanism. 
By means of equilibrium simulations the competition between strain and binding energy 
was found to yield regular stripe patterns, similar to the ones reported in experiments. 

But - as our lattice simulations proved - the edge diffusion barrier between regions 
of different particle types can already cause a stripe-like separation, without any con- 
sideration of an atomic mismatch. However, it seems clear that strain effects cause 
both, a restriction of the stripe width and a pronounced asymmetry in the behavior 
of the two different particle types. Such an asymmetry, where one particle type forms 
the backbone of a ramified island and the other adsorbate species gathers inbetween 
these branches is also reported in experimental studies (e.g. [37]). We have identified 
the different edge diffusion behavior of both adsorbate types as an important driving 
force for the observed island morphology. 

Despite the success in describing multi-component growth within our off-lattice 
model, several interesting questions still remain. One important issue is the influence 
of the lattice structure on the made observations. Especially the island morphology 
and the pattern formation in the equilibrium simulations are of interest here. It will be 
necessary to extend the method to more realistic lattice symmetries (e.g. fee). 

Also, the use of more material specific interaction potentials is an interesting task. 
Along with a realistic lattice structure this could allow for the prediction of phenomena 
in real systems. The application of so-called EAM (embedded atom method) many 
body potentials would be of particular relevance to metallic systems. Since most of 
the experimental systems exhibit misfits \e\ > 8%, dislocations are supposed to be of 



90 



some importance. The simulations should therefore be extended beyond the lattice 
based method. This would also facilitate investigations on the competition between 
misfit dislocations and surface alloying as possible strain relaxation mechanisms in the 
submonolayer regime. 
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Appendix A 

Pair-potentials used in this work 



In this work pair-potentials which depend on the distance \r^\ between two particles 
i and j of the crystal are used to compute the binding and transition energies for surface 
diffusion. 

The binding of two atoms - which is caused by changes in the electron density - is 
modeled by an approximation Uij of the interactions within the two-particle system. 
Many particle-interactions are neglected. This becomes problematic for the case of 
metallic systems where the bonding electrons are distributed over a large number of 
atoms and the binding-strength between two atoms is likely to depend highly on the 
local arrangement of the crystal (see e.g. [132]). 

However, pair-potentials are numerically easy to handle and provide the possibility 
of real-time calculations in KMC simulations. Since the aim of this work is not to model 
specific materials realistically but to gain general insight into relevant mechanisms of 
heteroepitaxial growth the advantages of pair-potentials outweigh their disadvantages. 



It has an attractive tail for large particle distance and reaches a minimum at the 
equilibrium distance 



where the depth of the potential is given by — Eij for the case n = 2m. The potential is 
zero at = and becomes strongly repulsive with decreasing particle distance [134]. 

The most common form of the Lennard- Jones potential is the 12, 6 form, where the 
attractive part oc 1/rf - is motivated by dipole-dipole interactions due to fluctuating 



A.1 Lennard-Jones potential 



The Lennard-Jones n, m potential [133] for two particles i and j is given by 




(A.l) 




(A.2) 
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Figure A.l: Lennard- Jones 12,6 and 
10, 5 potential for Oij = 1. 
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dipoles. These weak interactions dominate the bonding character of rare gases. The 
repulsive oc l/rjj term represents the repulsion between the electron clouds of the two 
atoms. The exponent n = 12 is chosen just for practical reasons since it allows a 
very cheap computation of the potential (see e.g. [134]). According to equation (A. 2) 
ro oc <Jij and the misfit e can be applied to the system by an appropriate choice of a if 
for substrate-substrate interaction we choose — a s — 1, the adsorbate-adsorbate 
interaction is given according to = a a = 1 + e and finally we set er^ = a as = l + e/2 
for the adsorbate-substrate interaction. 



Another pair-potential used in this work is the Morse potential [135] given by 



Like the Lennard- Jones potential it is attractive at large particle distances and 
reaches a minimum at the equilibrium distance r = a^, where the depth of the potential 
is given by —E^. 

In contrast to the Lennard- Jones potential the Morse potential remains finite in the 
repulsive part — > although it may become very large provided e aaiJ » 2. In this 
work the Morse potential is used for values of a between 4.0 and 7.0. As figure A. 2 
displays the potential becomes the steeper both for repulsive and attractive part, the 
larger the value of a is. In this way larger values of a lead to a better localization of 
the particles in the crystal. 

By using the experimental values for the energy of vaporization, the lattice constant 
and the compressibility, the Morse potential was fitted to a number of fee and bec 
metals [136]. This leads to values for a between 1.0 and 3.0. However, as mentioned 
above the application of pair-potentials to metallic systems is problematic. Therefore, 
fitted potentials give only a rough estimate of the binding energies in the bulk lattice. 



A.2 Morse potential 
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(A.3) 
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Figure A. 2: Morse potential for by = 
1 and different values of a. 



1 1.5 2 2.5 

ij 



In order to introduce the misfit e into an heteroepitaxial system we choose = 
a s — 1, (Tij — a a — 1 + e and = a as = 1 +e/2 for the substrate-substrate, adsorbate- 
adsorbate and adsorbate-substrate interaction, respectively. 

A.3 Simple cubic lattice 




Figure A.3: Contour plot of a cubic 
Lennard- Jones 12,6 potential with z = 
0.7. Darker regions correspond to lower 
values of U C vh- 
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Due to the isotropy of the central force potentials (A.l) and (A.3) particles arrange 
into triangular or fee lattices for the 2d and 3d case, respectively. However, sometimes a 
simple cubic lattice structure is favorable for growth simulations. Although this lattice 
structure does not represent any relevant material it has some advantages for basic 
examinations. Due to the lower coordination number less particles have to be taken 
into account for energy calculations than in a closed-packed lattice. There are only 
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four possible in-plane transition sites per particle, in comparison with six transition 
sites in case of, e.g., the fee lattice. This makes faster calculations of activation energies 
possible and allows for larger system sizes and more simulation runs. 

In order to stabilize a simple cubic lattice in our simulations an anisotropic potential, 



is used [11]. Uij denotes one of the central force potentials (A.l), (A. 3). The particle 
distance is given in Cartesian coordinates (x^-, y^-, Zij) where growth proceeds in 
z-direction. Figure A. 3 shows a contour plot for the Lennard- Jones 12, 6 type of the 
cubic potential (A.4) with z = 0.7. The binding sites around the central particle are 
arranged in a cubic symmetry with the local minima along the lattice directions. 




(A.4) 
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